Chapter 6 Diversity analysis

6.1 Alpha diversity

# Calculate Hill numbers
richness <- genome_counts_filt %>%
  column_to_rownames(var = "genome") %>%
  dplyr::select(where(~ !all(. == 0))) %>%
  hilldiv(., q = 0) %>%
  t() %>%
  as.data.frame() %>%
  dplyr::rename(richness = 1) %>%
  rownames_to_column(var = "sample")

neutral <- genome_counts_filt %>%
  column_to_rownames(var = "genome") %>%
  dplyr::select(where(~ !all(. == 0))) %>%
  hilldiv(., q = 1) %>%
  t() %>%
  as.data.frame() %>%
  dplyr::rename(neutral = 1) %>%
  rownames_to_column(var = "sample")

phylogenetic <- genome_counts_filt %>%
  column_to_rownames(var = "genome") %>%
  dplyr::select(where(~ !all(. == 0))) %>%
  hilldiv(., q = 1, tree = genome_tree) %>%
  t() %>%
  as.data.frame() %>%
  dplyr::rename(phylogenetic = 1) %>%
  rownames_to_column(var = "sample")

# Aggregate basal GIFT into elements
dist <- genome_gifts %>%
  to.elements(., GIFT_db3) %>%
  traits2dist(., method = "gower")

functional <- genome_counts_filt %>%
  column_to_rownames(var = "genome") %>%
  dplyr::select(where(~ !all(. == 0))) %>%
  hilldiv(., q = 1, dist = dist) %>%
  t() %>%
  as.data.frame() %>%
  dplyr::rename(functional = 1) %>%
  rownames_to_column(var = "sample") %>%
  mutate(functional = if_else(is.nan(functional), 1, functional))

# Merge all metrics
alpha_div <- richness %>%
  full_join(neutral, by = join_by(sample == sample)) %>%
  full_join(phylogenetic, by = join_by(sample == sample)) %>%
  full_join(functional, by = join_by(sample == sample))

6.1.1 Wild samples

alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point=="0_Wild") %>%
  mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
      ggplot(aes(y = value, x = Population, group=Population, color=Population, fill=Population)) +
      geom_boxplot(outlier.shape = NA) +
      geom_jitter(alpha=0.5) +
      scale_color_manual(name="Population",
          breaks=c("Cold_wet","Hot_dry"),
          labels=c("Cold","Hot"),
          values=c('#008080', "#d57d2c")) +
      scale_fill_manual(name="Population",
          breaks=c("Cold_wet","Hot_dry"),
          labels=c("Cold","Hot"),
          values=c('#00808050', "#d57d2c50")) +
      facet_wrap(. ~ metric,scales = "free") +
      coord_cartesian(xlim = c(1, NA)) +
      stat_compare_means(size=3, label.x=.58) +
      theme_classic() +
        theme(
    strip.background = element_blank(),
    panel.grid.minor.x = element_line(size = .1, color = "grey"),
    axis.title.x = element_blank(),
    axis.title.y = element_text(size=10),
    axis.text.x = element_text(angle = 45, hjust = 1),
    # Increase plot size
    plot.title = element_text(size = 10),
    axis.text = element_text(size = 8),
    axis.title = element_text(size = 8)
      ) +
  ylab("Alpha diversity")

6.1.2 Acclimation samples

alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point=="1_Acclimation") %>%
  mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
      ggplot(aes(y = value, x = Population, group=Population, color=Population, fill=Population)) +
      geom_boxplot(outlier.shape = NA) +
      geom_jitter(alpha=0.5) +
      scale_color_manual(name="Population",
          breaks=c("Cold_wet","Hot_dry"),
          labels=c("Cold","Hot"),
          values=c('#008080', "#d57d2c")) +
      scale_fill_manual(name="Population",
          breaks=c("Cold_wet","Hot_dry"),
          labels=c("Cold","Hot"),
          values=c('#00808050', "#d57d2c50")) +
      facet_wrap(. ~ metric,scales = "free") +
      coord_cartesian(xlim = c(1, NA)) +
      stat_compare_means(size=3, label.x=.58) +
      theme_classic() +
        theme(
    strip.background = element_blank(),
    panel.grid.minor.x = element_line(size = .1, color = "grey"),
    axis.title.x = element_blank(),
    axis.title.y = element_text(size=10),
    axis.text.x = element_text(angle = 45, hjust = 1),
    # Increase plot size
    plot.title = element_text(size = 10),
    axis.text = element_text(size = 8),
    axis.title = element_text(size = 8)
      ) +
  ylab("Alpha diversity")

6.1.3 Post-Transplant_1 samples

bxp1<-alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point=="5_Post-FMT1") %>%
  mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
      ggplot(aes(y = value, x = type, group=type, color=type, fill=type)) +
      geom_boxplot(outlier.shape = NA) +
      geom_jitter(alpha=0.5) +
      scale_color_manual(name="Type",
          breaks=c("Control","Hot_control", "Treatment"),
          labels=c("Cold-Cold","Hot-Hot", "Cold-Hot"),
          values=c("#4477AA","#d57d2c", "#76b183")) +
      scale_fill_manual(name="Type",
          breaks=c("Control","Hot_control", "Treatment"),
          labels=c("Cold-Cold","Hot-Hot", "Cold-Hot"),
          values=c("#4477AA50","#d57d2c50","#76b18350")) +
        scale_x_discrete(labels = c("Control" = "Cold-Cold", "Hot_control" = "Hot-Hot", "Treatment" = "Cold-Hot")) +
      facet_wrap(. ~ metric,scales = "free") +
      coord_cartesian(xlim = c(1, NA)) +
      stat_compare_means(size=3, label.x=.7) +
      theme_classic() +
        theme(
    strip.background = element_blank(),
    panel.grid.minor.x = element_line(size = .1, color = "grey"),
    axis.title.x = element_blank(),
    axis.title.y = element_text(size=10),
    axis.text.x = element_text(angle = 45, hjust = 1),
    # Increase plot size
    plot.title = element_text(size = 10),
    axis.text = element_text(size = 8),
    axis.title = element_text(size = 8)
      ) +
  ylab("Alpha diversity")


# Add significance bars with p-values
bxp1 + 
  geom_signif(comparisons = list(c("Control", "Hot_control"), c("Control", "Treatment"), c("Hot_control", "Treatment")),
              map_signif_level = TRUE)

6.1.4 Post-Transplant_2 samples

bxp1<-alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point=="6_Post-FMT2") %>%
  mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
      ggplot(aes(y = value, x = type, group=type, color=type, fill=type)) +
      geom_boxplot(outlier.shape = NA) +
      geom_jitter(alpha=0.5) +
      scale_color_manual(name="Type",
          breaks=c("Control","Hot_control", "Treatment"),
          labels=c("Cold-control","Warm-control", "Cold-intervention"),
          values=c("#4477AA","#d57d2c", "#76b183")) +
      scale_fill_manual(name="Type",
          breaks=c("Control","Hot_control", "Treatment"),
          labels=c("Cold-control","Warm-control", "Cold-intervention"),
          values=c("#4477AA50","#d57d2c50","#76b18350")) +
          scale_x_discrete(labels = c("Control" = "Cold-Cold", "Hot_control" = "Hot-Hot", "Treatment" = "Cold-Hot")) +
      facet_wrap(. ~ metric,scales = "free") +
      coord_cartesian(xlim = c(1, NA)) +
      stat_compare_means(size=3, label.x=.7) +
      theme_classic() +
        theme(
    strip.background = element_blank(),
    panel.grid.minor.x = element_line(size = .1, color = "grey"),
    axis.title.x = element_blank(),
    axis.title.y = element_text(size=10),
    axis.text.x = element_text(angle = 45, hjust = 1),
    # Increase plot size
    plot.title = element_text(size = 10),
    axis.text = element_text(size = 8),
    axis.title = element_text(size = 8)
      ) +
  ylab("Alpha diversity")

# Add significance bars with p-values
bxp1 + 
  geom_signif(comparisons = list(c("Control", "Hot_control"), c("Control", "Treatment"), c("Hot_control", "Treatment")),
              map_signif_level = TRUE)

6.2 Beta diversity

beta_q0n <- genome_counts_filt %>%
  column_to_rownames(., "genome") %>%
  hillpair(., q = 0)

beta_q1n <- genome_counts_filt %>%
  column_to_rownames(., "genome") %>%
  hillpair(., q = 1)

beta_q1p <- genome_counts_filt %>%
  column_to_rownames(., "genome") %>%
  hillpair(., q = 1, tree = genome_tree)

beta_q1f <- genome_counts_filt %>%
  column_to_rownames(., "genome") %>%
  hillpair(., q = 1, dist = dist)

6.3 Permanovas

6.3.0.1 Load required data

6.3.1 1. Are the wild populations similar?

6.3.1.1 Wild: cold vs warm adapted lizards

6.3.1.1.1 Number of samples used

[1] 27

6.3.1.2 Richness

richness <- as.matrix(beta_q0n$S)
richness <- as.dist(richness[rownames(richness) %in% samples_to_keep,
               colnames(richness) %in% samples_to_keep])
betadisper(richness, subset_meta$Population) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq  Mean Sq      F N.Perm Pr(>F)
Groups     1 0.000012 0.000012 0.0012    999  0.978
Residuals 25 0.257281 0.010291                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
         Cold_wet Hot_dry
Cold_wet             0.98
Hot_dry   0.97302        

adonis2(richness ~ Population,
        data = subset_meta %>% arrange(match(Tube_code,labels(richness))),
        permutations = 999) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	1	1.542719	0.2095041	6.625717	0.001
Residual	25	5.820951	0.7904959	NA	NA
Total	26	7.363669	1.0000000	NA	NA

6.3.1.3 Neutral

neutral <- as.matrix(beta_q1n$S)
neutral <- as.dist(neutral[rownames(neutral) %in% samples_to_keep,
               colnames(neutral) %in% samples_to_keep])
betadisper(neutral, subset_meta$Population) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     1 0.000048 0.0000476 0.0044    999  0.951
Residuals 25 0.270114 0.0108046                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
         Cold_wet Hot_dry
Cold_wet            0.954
Hot_dry   0.94763        

adonis2(neutral ~ Population,
        data = subset_meta %>% arrange(match(Tube_code,labels(neutral))),
        permutations = 999) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	1	1.918266	0.2608511	8.822682	0.001
Residual	25	5.435610	0.7391489	NA	NA
Total	26	7.353876	1.0000000	NA	NA

6.3.1.4 Phylogenetic

phylo <- as.matrix(beta_q1p$S)
phylo <- as.dist(phylo[rownames(phylo) %in% samples_to_keep,
               colnames(phylo) %in% samples_to_keep])
betadisper(phylo, subset_meta$Population) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)
Groups     1 0.03585 0.035847 2.4912    999  0.133
Residuals 25 0.35973 0.014389                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
         Cold_wet Hot_dry
Cold_wet            0.126
Hot_dry   0.12705        

adonis2(phylo ~ Population,
        data = subset_meta %>% arrange(match(Tube_code,labels(phylo))),
        permutations = 999) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	1	0.3218613	0.2162815	6.899207	0.001
Residual	25	1.1662981	0.7837185	NA	NA
Total	26	1.4881594	1.0000000	NA	NA

6.3.1.5 Functional

func <- as.matrix(beta_q1f$S)
func <- as.dist(func[rownames(func) %in% samples_to_keep,
               colnames(func) %in% samples_to_keep])
betadisper(func, subset_meta$Population) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq  Mean Sq     F N.Perm Pr(>F)
Groups     1 0.019387 0.019387 1.653    999  0.199
Residuals 25 0.293200 0.011728                    

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
         Cold_wet Hot_dry
Cold_wet            0.205
Hot_dry   0.21033        

adonis2(func ~ Population,
        data = subset_meta %>% arrange(match(Tube_code,labels(func))),
        permutations = 999) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	1	0.0831048	0.1680538	5.05002	0.067
Residual	25	0.4114083	0.8319462	NA	NA
Total	26	0.4945131	1.0000000	NA	NA

beta_q0n_nmds_wild <- richness %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE, trace=FALSE) %>%
                vegan::scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(subset_meta, by = join_by(sample == Tube_code))

beta_q1n_nmds_wild <- neutral %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE, trace=FALSE) %>%
                vegan::scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(subset_meta, by = join_by(sample == Tube_code))

beta_q1p_nmds_wild <- phylo %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                vegan::scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(subset_meta, by = join_by(sample == Tube_code))

beta_q1f_nmds_wild <- func %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                vegan::scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(subset_meta, by = join_by(sample == Tube_code))

6.3.2 2. Effect of acclimation

6.3.2.0.1 Number of samples used

[1] 26

6.3.2.1 Richness

richness_accli <- as.matrix(beta_q0n$S)
richness_accli <- as.dist(richness_accli[rownames(richness_accli) %in% samples_to_keep_accli,
               colnames(richness_accli) %in% samples_to_keep_accli])
betadisper(richness_accli, subset_meta_accli$Population) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq  Mean Sq      F N.Perm Pr(>F)   
Groups     1 0.093187 0.093187 11.812    999  0.002 **
Residuals 24 0.189340 0.007889                        
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
          Cold_wet Hot_dry
Cold_wet             0.002
Hot_dry  0.0021532        

adonis2(richness_accli ~ Population,
        data = subset_meta_accli %>% arrange(match(Tube_code,labels(richness_accli))),
        permutations = 999) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	1	1.639630	0.1929088	5.736415	0.001
Residual	24	6.859879	0.8070912	NA	NA
Total	25	8.499509	1.0000000	NA	NA

6.3.2.2 Neutral

neutral_accli <- as.matrix(beta_q1n$S)
neutral_accli <- as.dist(neutral_accli[rownames(neutral_accli) %in% samples_to_keep_accli,
               colnames(neutral_accli) %in% samples_to_keep_accli])
betadisper(neutral_accli, subset_meta_accli$Population) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)  
Groups     1 0.05603 0.056026 4.1918    999  0.049 *
Residuals 24 0.32077 0.013365                       
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
         Cold_wet Hot_dry
Cold_wet            0.051
Hot_dry  0.051717        

adonis2(neutral_accli ~ Population,
        data = subset_meta_accli %>% arrange(match(Tube_code,labels(neutral_accli))),
        permutations = 999) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	1	1.961192	0.2493171	7.970889	0.001
Residual	24	5.905063	0.7506829	NA	NA
Total	25	7.866255	1.0000000	NA	NA

6.3.2.3 Phylogenetic

phylo_accli <- as.matrix(beta_q1p$S)
phylo_accli <- as.dist(phylo_accli[rownames(phylo_accli) %in% samples_to_keep_accli,
               colnames(phylo_accli) %in% samples_to_keep_accli])
betadisper(phylo_accli, subset_meta_accli$Population) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)
Groups     1 0.03637 0.036365 2.3087    999  0.136
Residuals 24 0.37804 0.015752                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
         Cold_wet Hot_dry
Cold_wet            0.133
Hot_dry   0.14172        

adonis2(phylo_accli ~ Population,
        data = subset_meta_accli %>% arrange(match(Tube_code,labels(phylo_accli))),
        permutations = 999) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	1	0.2113395	0.1379515	3.84066	0.007
Residual	24	1.3206449	0.8620485	NA	NA
Total	25	1.5319844	1.0000000	NA	NA

6.3.2.4 Functional

func_accli <- as.matrix(beta_q1f$S)
func_accli <- as.dist(func_accli[rownames(func_accli) %in% samples_to_keep_accli,
               colnames(func_accli) %in% samples_to_keep_accli])
betadisper(func_accli, subset_meta_accli$Population) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)
Groups     1 0.00409 0.004087 0.1789    999  0.673
Residuals 24 0.54821 0.022842                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
         Cold_wet Hot_dry
Cold_wet             0.68
Hot_dry   0.67607        

adonis2(func_accli ~ Population,
        data = subset_meta_accli %>% arrange(match(Tube_code,labels(func_accli))),
        permutations = 999) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	1	0.00401769	0.005973179	0.1442177	0.636
Residual	24	0.66860401	0.994026821	NA	NA
Total	25	0.67262170	1.000000000	NA	NA

beta_q0n_nmds_accli <- richness_accli %>%
  metaMDS(.,trymax = 500, k=2, verbosity=FALSE, trace=FALSE) %>%
  vegan::scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(subset_meta_accli, by = join_by(sample == Tube_code))

beta_q1n_nmds_accli <- neutral_accli %>%
  metaMDS(.,trymax = 500, k=2, verbosity=FALSE, trace=FALSE) %>%
  vegan::scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(subset_meta_accli, by = join_by(sample == Tube_code))

beta_q1p_nmds_accli <- phylo_accli %>%
  metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
  vegan::scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(subset_meta_accli, by = join_by(sample == Tube_code))

beta_q1f_nmds_accli <- func_accli %>%
  metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
  vegan::scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(subset_meta_accli, by = join_by(sample == Tube_code))

6.3.3 3. Comparison between Wild and Acclimation

6.3.3.0.1 Number of samples used

[1] 53

6.3.3.0.2 Richness

richness_accli1 <- as.matrix(beta_q0n$S)
richness_accli1 <- as.dist(richness_accli1[rownames(richness_accli1) %in% samples_to_keep_accli1,
               colnames(richness_accli1) %in% samples_to_keep_accli1])
betadisper(richness_accli1, subset_meta_accli1$Population) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)  
Groups     1 0.03286 0.032865 2.9698    999  0.097 .
Residuals 51 0.56438 0.011066                       
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
         Cold_wet Hot_dry
Cold_wet            0.089
Hot_dry   0.09089        

adonis2(richness_accli1 ~ Population*time_point,
        data = subset_meta_accli1 %>% arrange(match(Tube_code,labels(richness_accli1))),
        permutations = 999,
        strata = subset_meta_accli1 %>% arrange(match(Tube_code,labels(richness_accli1))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	3	3.762933	0.2288365	4.846784	0.001
Residual	49	12.680830	0.7711635	NA	NA
Total	52	16.443763	1.0000000	NA	NA

#Arrange of metadata dataframe
subset_meta_accli1_arrange <- column_to_rownames(subset_meta_accli1, "Tube_code")
subset_meta_accli1_arrange<-subset_meta_accli1_arrange[labels(richness_accli1),]

pairwise<-pairwise.adonis(richness_accli1,subset_meta_accli1_arrange$Population_time, perm=999)
pairwise%>%
  tt()

pairs	Df	SumsOfSqs	F.Model	R2	p.value	p.adjusted	sig
Cold_wet.1_Acclimation vs Hot_dry.1_Acclimation	1	1.6396299	5.736415	0.19290877	0.001	0.006	*
Cold_wet.1_Acclimation vs Cold_wet.0_Wild	1	0.5273732	1.906862	0.05462715	0.008	0.048	.
Cold_wet.1_Acclimation vs Hot_dry.0_Wild	1	1.5558412	5.190707	0.17782052	0.001	0.006	*
Hot_dry.1_Acclimation vs Cold_wet.0_Wild	1	1.8259388	8.319131	0.24968031	0.001	0.006	*
Hot_dry.1_Acclimation vs Hot_dry.0_Wild	1	0.3856333	1.736034	0.09788177	0.008	0.048	.
Cold_wet.0_Wild vs Hot_dry.0_Wild	1	1.5427188	6.625717	0.20950408	0.001	0.006	*

6.3.3.0.3 Neutral

neutral_accli1 <- as.matrix(beta_q1n$S)
neutral_accli1 <- as.dist(neutral_accli1[rownames(neutral_accli1) %in% samples_to_keep_accli1,
               colnames(neutral_accli1) %in% samples_to_keep_accli1])
betadisper(neutral_accli1, subset_meta_accli1$Population) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)
Groups     1 0.01431 0.014310 1.1177    999  0.324
Residuals 51 0.65296 0.012803                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
         Cold_wet Hot_dry
Cold_wet            0.312
Hot_dry   0.29539        

adonis2(neutral_accli1 ~ Population*time_point,
        data = subset_meta_accli1 %>% arrange(match(Tube_code,labels(neutral_accli1))),
        permutations = 999,
        strata = subset_meta_accli1 %>% arrange(match(Tube_code,labels(neutral_accli1))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	3	4.735485	0.2945657	6.820252	0.001
Residual	49	11.340673	0.7054343	NA	NA
Total	52	16.076158	1.0000000	NA	NA

pairwise<-pairwise.adonis(neutral_accli1,subset_meta_accli1_arrange$Population_time, perm=999)
pairwise%>%
  tt()

pairs	Df	SumsOfSqs	F.Model	R2	p.value	p.adjusted	sig
Cold_wet.1_Acclimation vs Hot_dry.1_Acclimation	1	1.9611917	7.970889	0.24931708	0.001	0.006	*
Cold_wet.1_Acclimation vs Cold_wet.0_Wild	1	0.7466371	3.064914	0.08498327	0.001	0.006	*
Cold_wet.1_Acclimation vs Hot_dry.0_Wild	1	2.1421077	8.176866	0.25412251	0.001	0.006	*
Hot_dry.1_Acclimation vs Cold_wet.0_Wild	1	2.0371666	10.078295	0.28730857	0.001	0.006	*
Hot_dry.1_Acclimation vs Hot_dry.0_Wild	1	0.6314393	3.060027	0.16054681	0.001	0.006	*
Cold_wet.0_Wild vs Hot_dry.0_Wild	1	1.9182663	8.822682	0.26085105	0.001	0.006	*

6.3.3.0.4 Phylogenetic

phylo_accli1 <- as.matrix(beta_q1p$S)
phylo_accli1 <- as.dist(phylo_accli1[rownames(phylo_accli1) %in% samples_to_keep_accli1,
               colnames(phylo_accli1) %in% samples_to_keep_accli1])
betadisper(phylo_accli1, subset_meta_accli1$Population) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq   Mean Sq     F N.Perm Pr(>F)
Groups     1 0.00001 0.0000111 6e-04    999  0.981
Residuals 51 0.89017 0.0174543                    

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
         Cold_wet Hot_dry
Cold_wet            0.982
Hot_dry   0.97998        

adonis2(phylo_accli1 ~ Population*time_point,
        data = subset_meta_accli1 %>% arrange(match(Tube_code,labels(phylo_accli1))),
        permutations = 999,
        strata = subset_meta_accli1 %>% arrange(match(Tube_code,labels(phylo_accli1))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	3	0.7622566	0.2345983	5.006223	0.001
Residual	49	2.4869430	0.7654017	NA	NA
Total	52	3.2491997	1.0000000	NA	NA

pairwise<-pairwise.adonis(phylo_accli1,subset_meta_accli1_arrange$Population_time, perm=999)
pairwise%>%
  tt()

pairs	Df	SumsOfSqs	F.Model	R2	p.value	p.adjusted	sig
Cold_wet.1_Acclimation vs Hot_dry.1_Acclimation	1	0.2113395	3.840660	0.1379515	0.011	0.066
Cold_wet.1_Acclimation vs Cold_wet.0_Wild	1	0.2076830	4.216762	0.1133028	0.002	0.012	.
Cold_wet.1_Acclimation vs Hot_dry.0_Wild	1	0.3635093	4.799786	0.1666605	0.001	0.006	*
Hot_dry.1_Acclimation vs Cold_wet.0_Wild	1	0.2121479	7.924059	0.2406769	0.002	0.012	.
Hot_dry.1_Acclimation vs Hot_dry.0_Wild	1	0.2092433	3.885515	0.1953942	0.001	0.006	*
Cold_wet.0_Wild vs Hot_dry.0_Wild	1	0.3218613	6.899207	0.2162815	0.001	0.006	*

6.3.3.0.5 Functional

func_accli1 <- as.matrix(beta_q1f$S)
func_accli1 <- as.dist(func_accli1[rownames(func_accli1) %in% samples_to_keep_accli1,
               colnames(func_accli1) %in% samples_to_keep_accli1])
betadisper(func_accli1, subset_meta_accli1$Population) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)
Groups     1 0.01374 0.013738 0.7923    999  0.399
Residuals 51 0.88435 0.017340                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
         Cold_wet Hot_dry
Cold_wet              0.4
Hot_dry    0.3776        

adonis2(func_accli1 ~ Population*time_point,
        data = subset_meta_accli1 %>% arrange(match(Tube_code,labels(func_accli1))),
        permutations = 999,
        strata = subset_meta_accli1 %>% arrange(match(Tube_code,labels(func_accli1))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	3	0.1024301	0.08662587	1.54908	0.194
Residual	49	1.0800123	0.91337413	NA	NA
Total	52	1.1824424	1.00000000	NA	NA

pairwise<-pairwise.adonis(func_accli1,subset_meta_accli1_arrange$Population_time, perm=999)
pairwise%>%
  tt()

pairs	Df	SumsOfSqs	F.Model	R2	p.value	p.adjusted	sig
Cold_wet.1_Acclimation vs Hot_dry.1_Acclimation	1	0.0040176901	0.144217742	0.0059731793	0.634	1.000
Cold_wet.1_Acclimation vs Cold_wet.0_Wild	1	0.0002775823	0.009969841	0.0003020251	0.704	1.000
Cold_wet.1_Acclimation vs Hot_dry.0_Wild	1	0.0871206378	3.762920746	0.1355376396	0.079	0.474
Hot_dry.1_Acclimation vs Cold_wet.0_Wild	1	0.0025235151	0.120315257	0.0047895600	0.606	1.000
Hot_dry.1_Acclimation vs Hot_dry.0_Wild	1	0.0409732757	4.066330477	0.2026444487	0.092	0.552
Cold_wet.0_Wild vs Hot_dry.0_Wild	1	0.0831047968	5.050019518	0.1680537849	0.060	0.360

beta_richness_nmds_accli1 <- richness_accli1 %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(subset_meta_accli1, by = c("sample" = "Tube_code"))

beta_neutral_nmds_accli1 <- neutral_accli1 %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(subset_meta_accli1, by = c("sample" = "Tube_code"))

beta_phylo_nmds_accli1 <- phylo_accli1 %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(subset_meta_accli1, by = join_by(sample == Tube_code))

beta_func_nmds_accli1 <- func_accli1 %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(subset_meta_accli1, by = join_by(sample == Tube_code))

6.3.4 4. Effect of FMT on microbiota community

6.3.4.1 Comparison between Acclimation vs Post-FMT1

6.3.4.1.1 Number of samples used

[1] 52

6.3.4.1.2 Richness

richness_post3 <- as.matrix(beta_q0n$S)
richness_post3 <- as.dist(richness_post3[rownames(richness_post3) %in% samples_to_keep_post3,
               colnames(richness_post3) %in% samples_to_keep_post3])
betadisper(richness_post3, subset_meta_post3$Population) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)    
Groups     1 0.10843 0.108427 25.578    999  0.001 ***
Residuals 50 0.21195 0.004239                         
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
           Cold_wet Hot_dry
Cold_wet              0.001
Hot_dry  6.0928e-06        

adonis2(richness_post3 ~ Population*time_point,
        data = subset_meta_post3 %>% arrange(match(Tube_code,labels(richness_post3))),
        permutations = 999,
        strata = subset_meta_post3 %>% arrange(match(Tube_code,labels(richness_post3))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	3	3.472191	0.193401	3.836375	0.001
Residual	48	14.481131	0.806599	NA	NA
Total	51	17.953321	1.000000	NA	NA

#Arrange of metadata dataframe
subset_meta_post3_arrange <- column_to_rownames(subset_meta_post3, "Tube_code")
subset_meta_post3_arrange<-subset_meta_post3_arrange[labels(richness_post3),]

pairwise <- pairwise.adonis(richness_post3, subset_meta_post3_arrange$type_time, perm=999)
pairwise%>%
  tt()

pairs	Df	SumsOfSqs	F.Model	R2	p.value	p.adjusted	sig
Control.1_Acclimation vs Treatment.1_Acclimation	1	0.3657243	1.123239	0.06966584	0.235	1.000
Control.1_Acclimation vs Hot_control.1_Acclimation	1	1.2605860	4.943717	0.24788342	0.001	0.015	.
Control.1_Acclimation vs Control.5_Post-FMT1	1	0.6804988	2.100611	0.12283837	0.002	0.030	.
Control.1_Acclimation vs Treatment.5_Post-FMT1	1	0.8481258	2.673290	0.16033371	0.001	0.015	.
Control.1_Acclimation vs Hot_control.5_Post-FMT1	1	1.1180661	3.809275	0.20252111	0.001	0.015	.
Treatment.1_Acclimation vs Hot_control.1_Acclimation	1	1.3606630	5.087152	0.24124415	0.001	0.015	.
Treatment.1_Acclimation vs Control.5_Post-FMT1	1	0.7216200	2.172734	0.11956009	0.001	0.015	.
Treatment.1_Acclimation vs Treatment.5_Post-FMT1	1	0.9551308	2.926054	0.16322910	0.001	0.015	.
Treatment.1_Acclimation vs Hot_control.5_Post-FMT1	1	1.2263345	4.039487	0.20157637	0.001	0.015	.
Hot_control.1_Acclimation vs Control.5_Post-FMT1	1	1.4319792	5.384836	0.25180628	0.001	0.015	.
Hot_control.1_Acclimation vs Treatment.5_Post-FMT1	1	0.8172413	3.194690	0.17558364	0.001	0.015	.
Hot_control.1_Acclimation vs Hot_control.5_Post-FMT1	1	0.5796135	2.441615	0.13239702	0.001	0.015	.
Control.5_Post-FMT1 vs Treatment.5_Post-FMT1	1	0.5615418	1.729004	0.10335366	0.014	0.210
Control.5_Post-FMT1 vs Hot_control.5_Post-FMT1	1	0.8438429	2.793772	0.14865413	0.001	0.015	.
Treatment.5_Post-FMT1 vs Hot_control.5_Post-FMT1	1	0.3734921	1.268929	0.07799710	0.100	1.000

6.3.4.1.3 Neutral

neutral_post3 <- as.matrix(beta_q1n$S)
neutral_post3 <- as.dist(neutral_post3[rownames(neutral_post3) %in% samples_to_keep_post3,
               colnames(neutral_post3) %in% samples_to_keep_post3])
betadisper(neutral_post3, subset_meta_post3$Population) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)   
Groups     1 0.09089 0.090889 12.898    999  0.002 **
Residuals 50 0.35233 0.007047                        
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
           Cold_wet Hot_dry
Cold_wet              0.001
Hot_dry  0.00074926        

adonis2(neutral_post3 ~ Population*time_point,
        data = subset_meta_post3 %>% arrange(match(Tube_code,labels(neutral_post3))),
        permutations = 999,
        strata = subset_meta_post3 %>% arrange(match(Tube_code,labels(neutral_post3))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	3	4.443476	0.2628643	5.705637	0.001
Residual	48	12.460591	0.7371357	NA	NA
Total	51	16.904067	1.0000000	NA	NA

pairwise <- pairwise.adonis(neutral_post3, subset_meta_post3_arrange$type_time, perm=999)
pairwise%>%
  tt()

pairs	Df	SumsOfSqs	F.Model	R2	p.value	p.adjusted	sig
Control.1_Acclimation vs Treatment.1_Acclimation	1	0.2759858	0.9928976	0.06208366	0.457	1.000
Control.1_Acclimation vs Hot_control.1_Acclimation	1	1.4146464	6.5078325	0.30257965	0.001	0.015	.
Control.1_Acclimation vs Control.5_Post-FMT1	1	0.8153894	3.0970603	0.17113610	0.001	0.015	.
Control.1_Acclimation vs Treatment.5_Post-FMT1	1	1.1809241	4.4856470	0.24265567	0.001	0.015	.
Control.1_Acclimation vs Hot_control.5_Post-FMT1	1	1.4321524	5.7774260	0.27806264	0.001	0.015	.
Treatment.1_Acclimation vs Hot_control.1_Acclimation	1	1.6347704	6.8326887	0.29925029	0.001	0.015	.
Treatment.1_Acclimation vs Control.5_Post-FMT1	1	0.9517634	3.3715700	0.17404733	0.001	0.015	.
Treatment.1_Acclimation vs Treatment.5_Post-FMT1	1	1.3127773	4.6298256	0.23585668	0.001	0.015	.
Treatment.1_Acclimation vs Hot_control.5_Post-FMT1	1	1.6713369	6.2395460	0.28056085	0.001	0.015	.
Hot_control.1_Acclimation vs Control.5_Post-FMT1	1	1.5409781	6.8338056	0.29928456	0.001	0.015	.
Hot_control.1_Acclimation vs Treatment.5_Post-FMT1	1	0.9133614	4.0964534	0.21451383	0.001	0.015	.
Hot_control.1_Acclimation vs Hot_control.5_Post-FMT1	1	0.6954835	3.2951234	0.17077493	0.001	0.015	.
Control.5_Post-FMT1 vs Treatment.5_Post-FMT1	1	0.6051778	2.2508491	0.13047758	0.016	0.240
Control.5_Post-FMT1 vs Hot_control.5_Post-FMT1	1	1.0528902	4.1436369	0.20570451	0.001	0.015	.
Treatment.5_Post-FMT1 vs Hot_control.5_Post-FMT1	1	0.4150076	1.6372683	0.09840968	0.058	0.870

6.3.4.1.4 Phylogenetic

phylo_post3 <- as.matrix(beta_q1p$S)
phylo_post3 <- as.dist(phylo_post3[rownames(phylo_post3) %in% samples_to_keep_post3,
               colnames(phylo_post3) %in% samples_to_keep_post3])
betadisper(phylo_post3, subset_meta_post3$Population) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)
Groups     1 0.04088 0.040882 2.9254    999  0.102
Residuals 50 0.69874 0.013975                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
         Cold_wet Hot_dry
Cold_wet            0.093
Hot_dry  0.093394        

adonis2(phylo_post3 ~ Population*time_point,
        data = subset_meta_post3 %>% arrange(match(Tube_code,labels(phylo_post3))),
        permutations = 999,
        strata = subset_meta_post3 %>% arrange(match(Tube_code,labels(phylo_post3))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	3	0.6468577	0.2180259	4.461035	0.004
Residual	48	2.3200272	0.7819741	NA	NA
Total	51	2.9668849	1.0000000	NA	NA

pairwise <- pairwise.adonis(phylo_post3, subset_meta_post3_arrange$type_time, perm=999)
pairwise%>%
  tt()

pairs	Df	SumsOfSqs	F.Model	R2	p.value	p.adjusted	sig
Control.1_Acclimation vs Treatment.1_Acclimation	1	0.05690611	0.7893300	0.04999136	0.559	1.000
Control.1_Acclimation vs Hot_control.1_Acclimation	1	0.10866997	3.0547314	0.16919285	0.010	0.150
Control.1_Acclimation vs Control.5_Post-FMT1	1	0.11760287	2.9988032	0.16661126	0.018	0.270
Control.1_Acclimation vs Treatment.5_Post-FMT1	1	0.09624422	1.7247115	0.10968160	0.134	1.000
Control.1_Acclimation vs Hot_control.5_Post-FMT1	1	0.18586788	4.1904227	0.21836011	0.004	0.060
Treatment.1_Acclimation vs Hot_control.1_Acclimation	1	0.23108846	4.0521838	0.20208192	0.002	0.030	.
Treatment.1_Acclimation vs Control.5_Post-FMT1	1	0.26358465	4.3608960	0.21417997	0.005	0.075
Treatment.1_Acclimation vs Treatment.5_Post-FMT1	1	0.25319427	3.2738422	0.17915456	0.045	0.675
Treatment.1_Acclimation vs Hot_control.5_Post-FMT1	1	0.39050120	5.9837393	0.27218933	0.001	0.015	.
Hot_control.1_Acclimation vs Control.5_Post-FMT1	1	0.14203376	5.4200212	0.25303529	0.001	0.015	.
Hot_control.1_Acclimation vs Treatment.5_Post-FMT1	1	0.09666753	2.3682173	0.13635351	0.020	0.300
Hot_control.1_Acclimation vs Hot_control.5_Post-FMT1	1	0.09252600	2.9824958	0.15711821	0.006	0.090
Control.5_Post-FMT1 vs Treatment.5_Post-FMT1	1	0.01842535	0.4144162	0.02688498	0.774	1.000
Control.5_Post-FMT1 vs Hot_control.5_Post-FMT1	1	0.05987967	1.7387847	0.09802164	0.095	1.000
Treatment.5_Post-FMT1 vs Hot_control.5_Post-FMT1	1	0.03212966	0.6477782	0.04139746	0.708	1.000

6.3.4.1.5 Functional

func_post3 <- as.matrix(beta_q1f$S)
func_post3 <- as.dist(func_post3[rownames(func_post3) %in% samples_to_keep_post3,
               colnames(func_post3) %in% samples_to_keep_post3])
betadisper(func_post3, subset_meta_post3$Population) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     1 0.00116 0.0011552 0.0615    999  0.798
Residuals 50 0.93946 0.0187892                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
         Cold_wet Hot_dry
Cold_wet            0.791
Hot_dry   0.80518        

adonis2(func_post3 ~ Population*time_point,
        data = subset_meta_post3 %>% arrange(match(Tube_code,labels(func_post3))),
        permutations = 999,
        strata = subset_meta_post3 %>% arrange(match(Tube_code,labels(func_post3))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	3	0.02861503	0.01961684	0.3201497	0.424
Residual	48	1.43008228	0.98038316	NA	NA
Total	51	1.45869731	1.00000000	NA	NA

pairwise <- pairwise.adonis(func_post3, subset_meta_post3_arrange$type_time, perm=999)
pairwise%>%
  tt()

pairs	Df	SumsOfSqs	F.Model	R2	p.value	p.adjusted	sig
Control.1_Acclimation vs Treatment.1_Acclimation	1	0.022076664	0.65002192	0.041534889	0.400	1.000
Control.1_Acclimation vs Hot_control.1_Acclimation	1	0.016749605	0.52091146	0.033561912	0.475	1.000
Control.1_Acclimation vs Control.5_Post-FMT1	1	-0.008325555	-0.22800110	-0.015434681	0.902	1.000
Control.1_Acclimation vs Treatment.5_Post-FMT1	1	0.126527918	3.45550519	0.197960767	0.065	0.975
Control.1_Acclimation vs Hot_control.5_Post-FMT1	1	0.051008429	1.23951838	0.076327287	0.316	1.000
Treatment.1_Acclimation vs Hot_control.1_Acclimation	1	0.001022590	0.05430389	0.003382513	0.639	1.000
Treatment.1_Acclimation vs Control.5_Post-FMT1	1	0.002157067	0.09411569	0.005847832	0.619	1.000
Treatment.1_Acclimation vs Treatment.5_Post-FMT1	1	0.056602363	2.56037069	0.145803909	0.148	1.000
Treatment.1_Acclimation vs Hot_control.5_Post-FMT1	1	0.009569124	0.35095521	0.021463896	0.506	1.000
Hot_control.1_Acclimation vs Control.5_Post-FMT1	1	-0.001745663	-0.08225018	-0.005167199	0.703	1.000
Hot_control.1_Acclimation vs Treatment.5_Post-FMT1	1	0.057758674	2.84545622	0.159449901	0.157	1.000
Hot_control.1_Acclimation vs Hot_control.5_Post-FMT1	1	0.005575266	0.21803560	0.013444020	0.557	1.000
Control.5_Post-FMT1 vs Treatment.5_Post-FMT1	1	0.119540855	4.84764704	0.244242909	0.076	1.000
Control.5_Post-FMT1 vs Hot_control.5_Post-FMT1	1	0.052587837	1.77308932	0.099762584	0.212	1.000
Treatment.5_Post-FMT1 vs Hot_control.5_Post-FMT1	1	0.012980354	0.44307662	0.028690955	0.464	1.000

beta_richness_nmds_post3 <- richness_post3 %>%
  metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
  scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(subset_meta_post3, by = c("sample" = "Tube_code"))

beta_neutral_nmds_post3 <- neutral_post3 %>%
  metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
  scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(subset_meta_post3, by = c("sample" = "Tube_code"))

beta_phylo_nmds_post3 <- phylo_post3 %>%
  metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
  scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(subset_meta_post3, by = join_by(sample == Tube_code))

beta_func_nmds_post3 <- func_post3 %>%
  metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
  scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(subset_meta_post3, by = join_by(sample == Tube_code))

6.3.4.2 Comparison between Acclimation vs Post-FMT2

6.3.4.2.1 Number of samples used

[1] 53

6.3.4.2.2 Richness

richness_post4 <- as.matrix(beta_q0n$S)
richness_post4 <- as.dist(richness_post4[rownames(richness_post4) %in% samples_to_keep_post4,
               colnames(richness_post4) %in% samples_to_keep_post4])
betadisper(richness_post4, subset_meta_post4$Population) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)    
Groups     1 0.07832 0.078322 11.371    999  0.001 ***
Residuals 51 0.35129 0.006888                         
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
          Cold_wet Hot_dry
Cold_wet             0.001
Hot_dry  0.0014303        

adonis2(richness_post4 ~ Population*time_point,
        data = subset_meta_post4 %>% arrange(match(Tube_code,labels(richness_post4))),
        permutations = 999,
        strata = subset_meta_post4 %>% arrange(match(Tube_code,labels(richness_post4))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	3	3.279021	0.1937289	3.924535	0.001
Residual	49	13.646799	0.8062711	NA	NA
Total	52	16.925820	1.0000000	NA	NA

#Arrange of metadata dataframe
subset_meta_post4_arrange <- column_to_rownames(subset_meta_post4, "Tube_code")
subset_meta_post4_arrange<-subset_meta_post4_arrange[labels(richness_post4),]

pairwise <- pairwise.adonis(richness_post4, subset_meta_post4_arrange$type_time, perm=999)
pairwise%>%
  tt()

pairs	Df	SumsOfSqs	F.Model	R2	p.value	p.adjusted	sig
Control.1_Acclimation vs Treatment.1_Acclimation	1	0.3657243	1.123239	0.06966584	0.246	1.000
Control.1_Acclimation vs Hot_control.1_Acclimation	1	1.2605860	4.943717	0.24788342	0.001	0.015	.
Control.1_Acclimation vs Treatment.6_Post-FMT2	1	0.8072604	2.940901	0.16392161	0.001	0.015	.
Control.1_Acclimation vs Control.6_Post-FMT2	1	0.4817387	1.660775	0.09968176	0.028	0.420
Control.1_Acclimation vs Hot_control.6_Post-FMT2	1	1.1179704	3.885459	0.20573812	0.001	0.015	.
Treatment.1_Acclimation vs Hot_control.1_Acclimation	1	1.3606630	5.087152	0.24124415	0.001	0.015	.
Treatment.1_Acclimation vs Treatment.6_Post-FMT2	1	0.9130048	3.195028	0.16645080	0.001	0.015	.
Treatment.1_Acclimation vs Control.6_Post-FMT2	1	0.5959230	1.984036	0.11032208	0.001	0.015	.
Treatment.1_Acclimation vs Hot_control.6_Post-FMT2	1	1.2747787	4.275366	0.21086503	0.001	0.015	.
Hot_control.1_Acclimation vs Treatment.6_Post-FMT2	1	0.6397330	2.913695	0.15405213	0.001	0.015	.
Hot_control.1_Acclimation vs Control.6_Post-FMT2	1	1.4575447	6.224524	0.28007456	0.002	0.030	.
Hot_control.1_Acclimation vs Hot_control.6_Post-FMT2	1	0.3276169	1.412318	0.08111028	0.037	0.555
Treatment.6_Post-FMT2 vs Control.6_Post-FMT2	1	0.6463814	2.560441	0.13795154	0.001	0.015	.
Treatment.6_Post-FMT2 vs Hot_control.6_Post-FMT2	1	0.4796256	1.916520	0.10696943	0.001	0.015	.
Control.6_Post-FMT2 vs Hot_control.6_Post-FMT2	1	1.1305044	4.268317	0.21059061	0.001	0.015	.

6.3.4.2.3 Neutral

neutral_post4 <- as.matrix(beta_q1n$S)
neutral_post4 <- as.dist(neutral_post4[rownames(neutral_post4) %in% samples_to_keep_post4,
               colnames(neutral_post4) %in% samples_to_keep_post4])
betadisper(neutral_post4, subset_meta_post4$Population) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)    
Groups     1 0.07811 0.078108 9.6342    999  0.001 ***
Residuals 51 0.41347 0.008107                         
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
          Cold_wet Hot_dry
Cold_wet             0.002
Hot_dry  0.0031133        

adonis2(neutral_post4 ~ Population*time_point,
        data = subset_meta_post4 %>% arrange(match(Tube_code,labels(neutral_post4))),
        permutations = 999,
        strata = subset_meta_post4 %>% arrange(match(Tube_code,labels(neutral_post4))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	3	3.823089	0.2357666	5.038847	0.001
Residual	49	12.392477	0.7642334	NA	NA
Total	52	16.215567	1.0000000	NA	NA

pairwise <- pairwise.adonis(neutral_post4, subset_meta_post4_arrange$type_time, perm=999)
pairwise%>%
  tt()

pairs	Df	SumsOfSqs	F.Model	R2	p.value	p.adjusted	sig
Control.1_Acclimation vs Treatment.1_Acclimation	1	0.2759858	0.9928976	0.06208366	0.453	1.000
Control.1_Acclimation vs Hot_control.1_Acclimation	1	1.4146464	6.5078325	0.30257965	0.001	0.015	.
Control.1_Acclimation vs Treatment.6_Post-FMT2	1	1.1524353	4.9536068	0.24825621	0.001	0.015	.
Control.1_Acclimation vs Control.6_Post-FMT2	1	0.4748999	1.9609749	0.11561688	0.002	0.030	.
Control.1_Acclimation vs Hot_control.6_Post-FMT2	1	1.3503168	5.4081420	0.26499923	0.001	0.015	.
Treatment.1_Acclimation vs Hot_control.1_Acclimation	1	1.6347704	6.8326887	0.29925029	0.001	0.015	.
Treatment.1_Acclimation vs Treatment.6_Post-FMT2	1	1.3540292	5.3398081	0.25022756	0.001	0.015	.
Treatment.1_Acclimation vs Control.6_Post-FMT2	1	0.6311089	2.4041625	0.13063146	0.002	0.030	.
Treatment.1_Acclimation vs Hot_control.6_Post-FMT2	1	1.6125755	5.9825981	0.27215155	0.001	0.015	.
Hot_control.1_Acclimation vs Treatment.6_Post-FMT2	1	0.6202327	3.1519868	0.16457754	0.001	0.015	.
Hot_control.1_Acclimation vs Control.6_Post-FMT2	1	1.5701179	7.6327037	0.32297209	0.001	0.015	.
Hot_control.1_Acclimation vs Hot_control.6_Post-FMT2	1	0.3634438	1.7083388	0.09647087	0.035	0.525
Treatment.6_Post-FMT2 vs Control.6_Post-FMT2	1	1.0227481	4.6483346	0.22511910	0.001	0.015	.
Treatment.6_Post-FMT2 vs Hot_control.6_Post-FMT2	1	0.5010202	2.2065321	0.12119453	0.001	0.015	.
Control.6_Post-FMT2 vs Hot_control.6_Post-FMT2	1	1.3619424	5.7710313	0.26507845	0.001	0.015	.

6.3.4.2.4 Phylogenetic

phylo_post4 <- as.matrix(beta_q1p$S)
phylo_post4 <- as.dist(phylo_post4[rownames(phylo_post4) %in% samples_to_keep_post4,
               colnames(phylo_post4) %in% samples_to_keep_post4])
betadisper(phylo_post4, subset_meta_post4$Population) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq     F N.Perm Pr(>F)  
Groups     1 0.03098 0.030984 2.885    999  0.089 .
Residuals 51 0.54772 0.010740                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
         Cold_wet Hot_dry
Cold_wet             0.08
Hot_dry    0.0955        

adonis2(phylo_post4 ~ Population*time_point,
        data = subset_meta_post4 %>% arrange(match(Tube_code,labels(phylo_post4))),
        permutations = 999,
        strata = subset_meta_post4 %>% arrange(match(Tube_code,labels(phylo_post4))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	3	0.6466392	0.2442162	5.277784	0.001
Residual	49	2.0011759	0.7557838	NA	NA
Total	52	2.6478151	1.0000000	NA	NA

pairwise <- pairwise.adonis(phylo_post4, subset_meta_post4_arrange$type_time, perm=999)
pairwise%>%
  tt()

pairs	Df	SumsOfSqs	F.Model	R2	p.value	p.adjusted	sig
Control.1_Acclimation vs Treatment.1_Acclimation	1	0.05690611	0.789330	0.04999136	0.542	1.000
Control.1_Acclimation vs Hot_control.1_Acclimation	1	0.10866997	3.054731	0.16919285	0.010	0.150
Control.1_Acclimation vs Treatment.6_Post-FMT2	1	0.15591805	4.379209	0.22597458	0.006	0.090
Control.1_Acclimation vs Control.6_Post-FMT2	1	0.07099742	1.879364	0.11134092	0.111	1.000
Control.1_Acclimation vs Hot_control.6_Post-FMT2	1	0.18682367	4.878754	0.24542555	0.002	0.030	.
Treatment.1_Acclimation vs Hot_control.1_Acclimation	1	0.23108846	4.052184	0.20208192	0.007	0.105
Treatment.1_Acclimation vs Treatment.6_Post-FMT2	1	0.36496892	6.396667	0.28560797	0.001	0.015	.
Treatment.1_Acclimation vs Control.6_Post-FMT2	1	0.22628210	3.829222	0.19311005	0.024	0.360
Treatment.1_Acclimation vs Hot_control.6_Post-FMT2	1	0.34830814	5.846334	0.26761166	0.001	0.015	.
Hot_control.1_Acclimation vs Treatment.6_Post-FMT2	1	0.10002871	4.383624	0.21505615	0.001	0.015	.
Hot_control.1_Acclimation vs Control.6_Post-FMT2	1	0.12577510	5.060129	0.24027055	0.001	0.015	.
Hot_control.1_Acclimation vs Hot_control.6_Post-FMT2	1	0.06334378	2.499774	0.13512455	0.022	0.330
Treatment.6_Post-FMT2 vs Control.6_Post-FMT2	1	0.05927454	2.382025	0.12958449	0.021	0.315
Treatment.6_Post-FMT2 vs Hot_control.6_Post-FMT2	1	0.06906280	2.722460	0.14541146	0.005	0.075
Control.6_Post-FMT2 vs Hot_control.6_Post-FMT2	1	0.11081709	4.043656	0.20174244	0.001	0.015	.

6.3.4.2.5 Functional

func_post4 <- as.matrix(beta_q1f$S)
func_post4 <- as.dist(func_post4[rownames(func_post4) %in% samples_to_keep_post4,
               colnames(func_post4) %in% samples_to_keep_post4])
betadisper(func_post4, subset_meta_post4$Population) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq   Mean Sq     F N.Perm Pr(>F)
Groups     1 0.00006 0.0000601 0.003    999  0.966
Residuals 51 1.03808 0.0203544                    

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
         Cold_wet Hot_dry
Cold_wet            0.965
Hot_dry   0.95688        

adonis2(func_post4 ~ Population*time_point,
        data = subset_meta_post4 %>% arrange(match(Tube_code,labels(func_post4))),
        permutations = 999,
        strata = subset_meta_post4 %>% arrange(match(Tube_code,labels(func_post4))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	3	0.02585374	0.01884827	0.313769	0.537
Residual	49	1.34582345	0.98115173	NA	NA
Total	52	1.37167719	1.00000000	NA	NA

pairwise <- pairwise.adonis(func_post4, subset_meta_post4_arrange$type_time, perm=999)
pairwise%>%
  tt()

pairs	Df	SumsOfSqs	F.Model	R2	p.value	p.adjusted	sig
Control.1_Acclimation vs Treatment.1_Acclimation	1	0.0220766636	0.650021923	0.0415348890	0.439	1
Control.1_Acclimation vs Hot_control.1_Acclimation	1	0.0167496054	0.520911459	0.0335619117	0.457	1
Control.1_Acclimation vs Treatment.6_Post-FMT2	1	0.0376809680	1.119650309	0.0694587220	0.298	1
Control.1_Acclimation vs Control.6_Post-FMT2	1	0.0292200956	0.920511083	0.0578191917	0.386	1
Control.1_Acclimation vs Hot_control.6_Post-FMT2	1	0.0133433458	0.270473784	0.0177122064	0.496	1
Treatment.1_Acclimation vs Hot_control.1_Acclimation	1	0.0010225901	0.054303886	0.0033825127	0.676	1
Treatment.1_Acclimation vs Treatment.6_Post-FMT2	1	-0.0005177706	-0.025585400	-0.0016016487	0.724	1
Treatment.1_Acclimation vs Control.6_Post-FMT2	1	0.0013301207	0.072110871	0.0044867082	0.610	1
Treatment.1_Acclimation vs Hot_control.6_Post-FMT2	1	0.0060959077	0.174487757	0.0107878382	0.576	1
Hot_control.1_Acclimation vs Treatment.6_Post-FMT2	1	0.0010345754	0.055797964	0.0034752533	0.629	1
Hot_control.1_Acclimation vs Control.6_Post-FMT2	1	-0.0001056284	-0.006306177	-0.0003942915	0.706	1
Hot_control.1_Acclimation vs Hot_control.6_Post-FMT2	1	0.0017235602	0.051851181	0.0032302306	0.773	1
Treatment.6_Post-FMT2 vs Control.6_Post-FMT2	1	-0.0080428882	-0.442986255	-0.0284750185	0.848	1
Treatment.6_Post-FMT2 vs Hot_control.6_Post-FMT2	1	-0.0011796256	-0.034047378	-0.0021324990	0.889	1
Control.6_Post-FMT2 vs Hot_control.6_Post-FMT2	1	0.0036300838	0.110487573	0.0068581148	0.706	1

beta_richness_nmds_post4 <- richness_post4 %>%
  metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
  scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(subset_meta_post4, by = c("sample" = "Tube_code"))

beta_neutral_nmds_post4 <- neutral_post4 %>%
  metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
  scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(subset_meta_post4, by = c("sample" = "Tube_code"))

beta_phylo_nmds_post4 <- phylo_post4 %>%
  metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
  scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(subset_meta_post4, by = join_by(sample == Tube_code))

beta_func_nmds_post4 <- func_post4 %>%
  metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
  scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(subset_meta_post4, by = join_by(sample == Tube_code))

6.3.4.3 Comparison between Acclimation vs Post-FMT1 and Post-FMT2

6.3.4.3.1 Number of samples used

[1] 79

6.3.4.3.2 Richness

richness_post6 <- as.matrix(beta_q0n$S)
richness_post6 <- as.dist(richness_post6[rownames(richness_post6) %in% samples_to_keep_post6,
               colnames(richness_post6) %in% samples_to_keep_post6])
betadisper(richness_post6, subset_meta_post6$Population) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)   
Groups     1 0.09792 0.097923 18.607    999  0.002 **
Residuals 77 0.40523 0.005263                        
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
          Cold_wet Hot_dry
Cold_wet             0.001
Hot_dry  4.714e-05        

adonis2(richness_post6 ~ type*time_point,
        data = subset_meta_post6 %>% arrange(match(Tube_code,labels(richness_post6))),
        permutations = 999,
        strata = subset_meta_post6 %>% arrange(match(Tube_code,labels(richness_post6))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	8	6.395467	0.2451322	2.841433	0.001
Residual	70	19.694403	0.7548678	NA	NA
Total	78	26.089870	1.0000000	NA	NA

#Arrange of metadata dataframe
subset_meta_post6_arrange <- column_to_rownames(subset_meta_post6, "Tube_code")
subset_meta_post6_arrange<-subset_meta_post6_arrange[labels(richness_post6),]

pairwise <- pairwise.adonis(richness_post6, subset_meta_post6_arrange$type_time, perm=999)
pairwise%>%
  tt()

pairs	Df	SumsOfSqs	F.Model	R2	p.value	p.adjusted	sig
Control.1_Acclimation vs Treatment.1_Acclimation	1	0.3657243	1.123239	0.06966584	0.248	1.000
Control.1_Acclimation vs Hot_control.1_Acclimation	1	1.2605860	4.943717	0.24788342	0.001	0.036	.
Control.1_Acclimation vs Control.5_Post-FMT1	1	0.6804988	2.100611	0.12283837	0.001	0.036	.
Control.1_Acclimation vs Treatment.5_Post-FMT1	1	0.8481258	2.673290	0.16033371	0.001	0.036	.
Control.1_Acclimation vs Hot_control.5_Post-FMT1	1	1.1180661	3.809275	0.20252111	0.001	0.036	.
Control.1_Acclimation vs Treatment.6_Post-FMT2	1	0.8072604	2.940901	0.16392161	0.001	0.036	.
Control.1_Acclimation vs Control.6_Post-FMT2	1	0.4817387	1.660775	0.09968176	0.026	0.936
Control.1_Acclimation vs Hot_control.6_Post-FMT2	1	1.1179704	3.885459	0.20573812	0.001	0.036	.
Treatment.1_Acclimation vs Hot_control.1_Acclimation	1	1.3606630	5.087152	0.24124415	0.001	0.036	.
Treatment.1_Acclimation vs Control.5_Post-FMT1	1	0.7216200	2.172734	0.11956009	0.001	0.036	.
Treatment.1_Acclimation vs Treatment.5_Post-FMT1	1	0.9551308	2.926054	0.16322910	0.001	0.036	.
Treatment.1_Acclimation vs Hot_control.5_Post-FMT1	1	1.2263345	4.039487	0.20157637	0.001	0.036	.
Treatment.1_Acclimation vs Treatment.6_Post-FMT2	1	0.9130048	3.195028	0.16645080	0.001	0.036	.
Treatment.1_Acclimation vs Control.6_Post-FMT2	1	0.5959230	1.984036	0.11032208	0.001	0.036	.
Treatment.1_Acclimation vs Hot_control.6_Post-FMT2	1	1.2747787	4.275366	0.21086503	0.001	0.036	.
Hot_control.1_Acclimation vs Control.5_Post-FMT1	1	1.4319792	5.384836	0.25180628	0.002	0.072
Hot_control.1_Acclimation vs Treatment.5_Post-FMT1	1	0.8172413	3.194690	0.17558364	0.001	0.036	.
Hot_control.1_Acclimation vs Hot_control.5_Post-FMT1	1	0.5796135	2.441615	0.13239702	0.001	0.036	.
Hot_control.1_Acclimation vs Treatment.6_Post-FMT2	1	0.6397330	2.913695	0.15405213	0.001	0.036	.
Hot_control.1_Acclimation vs Control.6_Post-FMT2	1	1.4575447	6.224524	0.28007456	0.001	0.036	.
Hot_control.1_Acclimation vs Hot_control.6_Post-FMT2	1	0.3276169	1.412318	0.08111028	0.038	1.000
Control.5_Post-FMT1 vs Treatment.5_Post-FMT1	1	0.5615418	1.729004	0.10335366	0.010	0.360
Control.5_Post-FMT1 vs Hot_control.5_Post-FMT1	1	0.8438429	2.793772	0.14865413	0.001	0.036	.
Control.5_Post-FMT1 vs Treatment.6_Post-FMT2	1	0.7628135	2.683925	0.14364890	0.002	0.072
Control.5_Post-FMT1 vs Control.6_Post-FMT2	1	0.3432605	1.148733	0.06698647	0.236	1.000
Control.5_Post-FMT1 vs Hot_control.6_Post-FMT2	1	1.1269580	3.799256	0.19188884	0.001	0.036	.
Treatment.5_Post-FMT1 vs Hot_control.5_Post-FMT1	1	0.3734921	1.268929	0.07799710	0.122	1.000
Treatment.5_Post-FMT1 vs Treatment.6_Post-FMT2	1	0.3571397	1.297184	0.07959561	0.139	1.000
Treatment.5_Post-FMT1 vs Control.6_Post-FMT2	1	0.7769467	2.670898	0.15114670	0.001	0.036	.
Treatment.5_Post-FMT1 vs Hot_control.6_Post-FMT2	1	0.6502360	2.253407	0.13060650	0.002	0.072
Hot_control.5_Post-FMT1 vs Treatment.6_Post-FMT2	1	0.4132091	1.616138	0.09174188	0.013	0.468
Hot_control.5_Post-FMT1 vs Control.6_Post-FMT2	1	1.0163992	3.760571	0.19030682	0.001	0.036	.
Hot_control.5_Post-FMT1 vs Hot_control.6_Post-FMT2	1	0.2732563	1.019281	0.05988979	0.416	1.000
Treatment.6_Post-FMT2 vs Control.6_Post-FMT2	1	0.6463814	2.560441	0.13795154	0.001	0.036	.
Treatment.6_Post-FMT2 vs Hot_control.6_Post-FMT2	1	0.4796256	1.916520	0.10696943	0.001	0.036	.
Control.6_Post-FMT2 vs Hot_control.6_Post-FMT2	1	1.1305044	4.268317	0.21059061	0.001	0.036	.

6.3.4.3.3 Neutral

neutral_post6 <- as.matrix(beta_q1n$S)
neutral_post6 <- as.dist(neutral_post6[rownames(neutral_post6) %in% samples_to_keep_post6,
               colnames(neutral_post6) %in% samples_to_keep_post6])
betadisper(neutral_post6, subset_meta_post6$Population) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)    
Groups     1 0.09524 0.095237 15.396    999  0.001 ***
Residuals 77 0.47632 0.006186                         
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
           Cold_wet Hot_dry
Cold_wet              0.001
Hot_dry  0.00018818        

adonis2(neutral_post6 ~ type*time_point,
        data = subset_meta_post6 %>% arrange(match(Tube_code,labels(neutral_post6))),
        permutations = 999,
        strata = subset_meta_post6 %>% arrange(match(Tube_code,labels(neutral_post6))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	8	7.736408	0.3121974	3.971673	0.001
Residual	70	17.044094	0.6878026	NA	NA
Total	78	24.780501	1.0000000	NA	NA

pairwise <- pairwise.adonis(neutral_post6, subset_meta_post6_arrange$type_time, perm=999)
pairwise%>%
  tt()

pairs	Df	SumsOfSqs	F.Model	R2	p.value	p.adjusted	sig
Control.1_Acclimation vs Treatment.1_Acclimation	1	0.2759858	0.9928976	0.06208366	0.460	1.000
Control.1_Acclimation vs Hot_control.1_Acclimation	1	1.4146464	6.5078325	0.30257965	0.001	0.036	.
Control.1_Acclimation vs Control.5_Post-FMT1	1	0.8153894	3.0970603	0.17113610	0.001	0.036	.
Control.1_Acclimation vs Treatment.5_Post-FMT1	1	1.1809241	4.4856470	0.24265567	0.001	0.036	.
Control.1_Acclimation vs Hot_control.5_Post-FMT1	1	1.4321524	5.7774260	0.27806264	0.001	0.036	.
Control.1_Acclimation vs Treatment.6_Post-FMT2	1	1.1524353	4.9536068	0.24825621	0.001	0.036	.
Control.1_Acclimation vs Control.6_Post-FMT2	1	0.4748999	1.9609749	0.11561688	0.006	0.216
Control.1_Acclimation vs Hot_control.6_Post-FMT2	1	1.3503168	5.4081420	0.26499923	0.001	0.036	.
Treatment.1_Acclimation vs Hot_control.1_Acclimation	1	1.6347704	6.8326887	0.29925029	0.002	0.072
Treatment.1_Acclimation vs Control.5_Post-FMT1	1	0.9517634	3.3715700	0.17404733	0.001	0.036	.
Treatment.1_Acclimation vs Treatment.5_Post-FMT1	1	1.3127773	4.6298256	0.23585668	0.001	0.036	.
Treatment.1_Acclimation vs Hot_control.5_Post-FMT1	1	1.6713369	6.2395460	0.28056085	0.001	0.036	.
Treatment.1_Acclimation vs Treatment.6_Post-FMT2	1	1.3540292	5.3398081	0.25022756	0.001	0.036	.
Treatment.1_Acclimation vs Control.6_Post-FMT2	1	0.6311089	2.4041625	0.13063146	0.001	0.036	.
Treatment.1_Acclimation vs Hot_control.6_Post-FMT2	1	1.6125755	5.9825981	0.27215155	0.001	0.036	.
Hot_control.1_Acclimation vs Control.5_Post-FMT1	1	1.5409781	6.8338056	0.29928456	0.001	0.036	.
Hot_control.1_Acclimation vs Treatment.5_Post-FMT1	1	0.9133614	4.0964534	0.21451383	0.001	0.036	.
Hot_control.1_Acclimation vs Hot_control.5_Post-FMT1	1	0.6954835	3.2951234	0.17077493	0.001	0.036	.
Hot_control.1_Acclimation vs Treatment.6_Post-FMT2	1	0.6202327	3.1519868	0.16457754	0.001	0.036	.
Hot_control.1_Acclimation vs Control.6_Post-FMT2	1	1.5701179	7.6327037	0.32297209	0.001	0.036	.
Hot_control.1_Acclimation vs Hot_control.6_Post-FMT2	1	0.3634438	1.7083388	0.09647087	0.035	1.000
Control.5_Post-FMT1 vs Treatment.5_Post-FMT1	1	0.6051778	2.2508491	0.13047758	0.013	0.468
Control.5_Post-FMT1 vs Hot_control.5_Post-FMT1	1	1.0528902	4.1436369	0.20570451	0.001	0.036	.
Control.5_Post-FMT1 vs Treatment.6_Post-FMT2	1	0.8908158	3.7146920	0.18842252	0.001	0.036	.
Control.5_Post-FMT1 vs Control.6_Post-FMT2	1	0.3860927	1.5521758	0.08843210	0.081	1.000
Control.5_Post-FMT1 vs Hot_control.6_Post-FMT2	1	1.3122237	5.1302726	0.24279254	0.001	0.036	.
Treatment.5_Post-FMT1 vs Hot_control.5_Post-FMT1	1	0.4150076	1.6372683	0.09840968	0.061	1.000
Treatment.5_Post-FMT1 vs Treatment.6_Post-FMT2	1	0.3157079	1.3252026	0.08117526	0.157	1.000
Treatment.5_Post-FMT1 vs Control.6_Post-FMT2	1	1.0579520	4.2700097	0.22158835	0.001	0.036	.
Treatment.5_Post-FMT1 vs Hot_control.6_Post-FMT2	1	0.7454015	2.9200493	0.16294873	0.002	0.072
Hot_control.5_Post-FMT1 vs Treatment.6_Post-FMT2	1	0.4377161	1.9421261	0.10824392	0.003	0.108
Hot_control.5_Post-FMT1 vs Control.6_Post-FMT2	1	1.3766597	5.8752789	0.26858075	0.001	0.036	.
Hot_control.5_Post-FMT1 vs Hot_control.6_Post-FMT2	1	0.3176516	1.3161367	0.07600637	0.184	1.000
Treatment.6_Post-FMT2 vs Control.6_Post-FMT2	1	1.0227481	4.6483346	0.22511910	0.001	0.036	.
Treatment.6_Post-FMT2 vs Hot_control.6_Post-FMT2	1	0.5010202	2.2065321	0.12119453	0.002	0.072
Control.6_Post-FMT2 vs Hot_control.6_Post-FMT2	1	1.3619424	5.7710313	0.26507845	0.001	0.036	.

6.3.4.3.4 Phylogenetic

phylo_post6 <- as.matrix(beta_q1p$S)
phylo_post6 <- as.dist(phylo_post6[rownames(phylo_post6) %in% samples_to_keep_post6,
               colnames(phylo_post6) %in% samples_to_keep_post6])
betadisper(phylo_post6, subset_meta_post6$Population) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)
Groups     1 0.02554 0.025541 2.1633    999  0.153
Residuals 77 0.90909 0.011806                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
         Cold_wet Hot_dry
Cold_wet            0.153
Hot_dry   0.14542        

adonis2(phylo_post6 ~ type*time_point,
        data = subset_meta_post6 %>% arrange(match(Tube_code,labels(phylo_post6))),
        permutations = 999,
        strata = subset_meta_post6 %>% arrange(match(Tube_code,labels(phylo_post6))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	8	1.055866	0.2692287	3.223651	0.001
Residual	70	2.865952	0.7307713	NA	NA
Total	78	3.921818	1.0000000	NA	NA

pairwise <- pairwise.adonis(phylo_post6, subset_meta_post6_arrange$type_time, perm=999)
pairwise%>%
  tt()

pairs	Df	SumsOfSqs	F.Model	R2	p.value	p.adjusted	sig
Control.1_Acclimation vs Treatment.1_Acclimation	1	0.05690611	0.7893300	0.04999136	0.506	1.000
Control.1_Acclimation vs Hot_control.1_Acclimation	1	0.10866997	3.0547314	0.16919285	0.012	0.432
Control.1_Acclimation vs Control.5_Post-FMT1	1	0.11760287	2.9988032	0.16661126	0.022	0.792
Control.1_Acclimation vs Treatment.5_Post-FMT1	1	0.09624422	1.7247115	0.10968160	0.161	1.000
Control.1_Acclimation vs Hot_control.5_Post-FMT1	1	0.18586788	4.1904227	0.21836011	0.006	0.216
Control.1_Acclimation vs Treatment.6_Post-FMT2	1	0.15591805	4.3792085	0.22597458	0.002	0.072
Control.1_Acclimation vs Control.6_Post-FMT2	1	0.07099742	1.8793639	0.11134092	0.101	1.000
Control.1_Acclimation vs Hot_control.6_Post-FMT2	1	0.18682367	4.8787543	0.24542555	0.001	0.036	.
Treatment.1_Acclimation vs Hot_control.1_Acclimation	1	0.23108846	4.0521838	0.20208192	0.005	0.180
Treatment.1_Acclimation vs Control.5_Post-FMT1	1	0.26358465	4.3608960	0.21417997	0.009	0.324
Treatment.1_Acclimation vs Treatment.5_Post-FMT1	1	0.25319427	3.2738422	0.17915456	0.041	1.000
Treatment.1_Acclimation vs Hot_control.5_Post-FMT1	1	0.39050120	5.9837393	0.27218933	0.002	0.072
Treatment.1_Acclimation vs Treatment.6_Post-FMT2	1	0.36496892	6.3966666	0.28560797	0.002	0.072
Treatment.1_Acclimation vs Control.6_Post-FMT2	1	0.22628210	3.8292220	0.19311005	0.019	0.684
Treatment.1_Acclimation vs Hot_control.6_Post-FMT2	1	0.34830814	5.8463335	0.26761166	0.001	0.036	.
Hot_control.1_Acclimation vs Control.5_Post-FMT1	1	0.14203376	5.4200212	0.25303529	0.001	0.036	.
Hot_control.1_Acclimation vs Treatment.5_Post-FMT1	1	0.09666753	2.3682173	0.13635351	0.019	0.684
Hot_control.1_Acclimation vs Hot_control.5_Post-FMT1	1	0.09252600	2.9824958	0.15711821	0.009	0.324
Hot_control.1_Acclimation vs Treatment.6_Post-FMT2	1	0.10002871	4.3836237	0.21505615	0.001	0.036	.
Hot_control.1_Acclimation vs Control.6_Post-FMT2	1	0.12577510	5.0601287	0.24027055	0.001	0.036	.
Hot_control.1_Acclimation vs Hot_control.6_Post-FMT2	1	0.06334378	2.4997737	0.13512455	0.015	0.540
Control.5_Post-FMT1 vs Treatment.5_Post-FMT1	1	0.01842535	0.4144162	0.02688498	0.791	1.000
Control.5_Post-FMT1 vs Hot_control.5_Post-FMT1	1	0.05987967	1.7387847	0.09802164	0.107	1.000
Control.5_Post-FMT1 vs Treatment.6_Post-FMT2	1	0.07917244	3.0180046	0.15869197	0.010	0.360
Control.5_Post-FMT1 vs Control.6_Post-FMT2	1	0.04335491	1.5335604	0.08746429	0.179	1.000
Control.5_Post-FMT1 vs Hot_control.6_Post-FMT2	1	0.10783045	3.7500438	0.18987521	0.001	0.036	.
Treatment.5_Post-FMT1 vs Hot_control.5_Post-FMT1	1	0.03212966	0.6477782	0.04139746	0.707	1.000
Treatment.5_Post-FMT1 vs Treatment.6_Post-FMT2	1	0.06393539	1.5651817	0.09448624	0.147	1.000
Treatment.5_Post-FMT1 vs Control.6_Post-FMT2	1	0.05265949	1.2240203	0.07544494	0.293	1.000
Treatment.5_Post-FMT1 vs Hot_control.6_Post-FMT2	1	0.09753501	2.2402429	0.12994265	0.012	0.432
Hot_control.5_Post-FMT1 vs Treatment.6_Post-FMT2	1	0.07228545	2.3279593	0.12701683	0.036	1.000
Hot_control.5_Post-FMT1 vs Control.6_Post-FMT2	1	0.11759094	3.5538444	0.18174658	0.001	0.036	.
Hot_control.5_Post-FMT1 vs Hot_control.6_Post-FMT2	1	0.06667255	1.9859527	0.11041687	0.108	1.000
Treatment.6_Post-FMT2 vs Control.6_Post-FMT2	1	0.05927454	2.3820253	0.12958449	0.026	0.936
Treatment.6_Post-FMT2 vs Hot_control.6_Post-FMT2	1	0.06906280	2.7224602	0.14541146	0.005	0.180
Control.6_Post-FMT2 vs Hot_control.6_Post-FMT2	1	0.11081709	4.0436561	0.20174244	0.002	0.072

6.3.4.3.5 Functional

func_post6 <- as.matrix(beta_q1f$S)
func_post6 <- as.dist(func_post6[rownames(func_post6) %in% samples_to_keep_post6,
               colnames(func_post6) %in% samples_to_keep_post6])
betadisper(func_post6, subset_meta_post6$Population) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)
Groups     1 0.00018 0.000175 0.0093    999  0.926
Residuals 77 1.45162 0.018852                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
         Cold_wet Hot_dry
Cold_wet            0.921
Hot_dry   0.92349        

adonis2(func_post6 ~ type*time_point,
        data = subset_meta_post6 %>% arrange(match(Tube_code,labels(func_post6))),
        permutations = 999,
        strata = subset_meta_post6 %>% arrange(match(Tube_code,labels(func_post6))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	8	0.1632318	0.07638509	0.7236452	0.56
Residual	70	1.9737271	0.92361491	NA	NA
Total	78	2.1369589	1.00000000	NA	NA

pairwise <- pairwise.adonis(func_post6, subset_meta_post6_arrange$type_time, perm=999)
pairwise%>%
  tt()

pairs	Df	SumsOfSqs	F.Model	R2	p.value	p.adjusted	sig
Control.1_Acclimation vs Treatment.1_Acclimation	1	0.0220766636	0.650021923	0.0415348890	0.400	1
Control.1_Acclimation vs Hot_control.1_Acclimation	1	0.0167496054	0.520911459	0.0335619117	0.448	1
Control.1_Acclimation vs Control.5_Post-FMT1	1	-0.0083255553	-0.228001097	-0.0154346814	0.888	1
Control.1_Acclimation vs Treatment.5_Post-FMT1	1	0.1265279179	3.455505193	0.1979607668	0.082	1
Control.1_Acclimation vs Hot_control.5_Post-FMT1	1	0.0510084290	1.239518381	0.0763272870	0.286	1
Control.1_Acclimation vs Treatment.6_Post-FMT2	1	0.0376809680	1.119650309	0.0694587220	0.305	1
Control.1_Acclimation vs Control.6_Post-FMT2	1	0.0292200956	0.920511083	0.0578191917	0.379	1
Control.1_Acclimation vs Hot_control.6_Post-FMT2	1	0.0133433458	0.270473784	0.0177122064	0.507	1
Treatment.1_Acclimation vs Hot_control.1_Acclimation	1	0.0010225901	0.054303886	0.0033825127	0.651	1
Treatment.1_Acclimation vs Control.5_Post-FMT1	1	0.0021570675	0.094115687	0.0058478321	0.635	1
Treatment.1_Acclimation vs Treatment.5_Post-FMT1	1	0.0566023632	2.560370692	0.1458039091	0.185	1
Treatment.1_Acclimation vs Hot_control.5_Post-FMT1	1	0.0095691236	0.350955208	0.0214638964	0.495	1
Treatment.1_Acclimation vs Treatment.6_Post-FMT2	1	-0.0005177706	-0.025585400	-0.0016016487	0.750	1
Treatment.1_Acclimation vs Control.6_Post-FMT2	1	0.0013301207	0.072110871	0.0044867082	0.606	1
Treatment.1_Acclimation vs Hot_control.6_Post-FMT2	1	0.0060959077	0.174487757	0.0107878382	0.576	1
Hot_control.1_Acclimation vs Control.5_Post-FMT1	1	-0.0017456629	-0.082250179	-0.0051671989	0.703	1
Hot_control.1_Acclimation vs Treatment.5_Post-FMT1	1	0.0577586745	2.845456220	0.1594499006	0.172	1
Hot_control.1_Acclimation vs Hot_control.5_Post-FMT1	1	0.0055752661	0.218035597	0.0134440201	0.530	1
Hot_control.1_Acclimation vs Treatment.6_Post-FMT2	1	0.0010345754	0.055797964	0.0034752533	0.671	1
Hot_control.1_Acclimation vs Control.6_Post-FMT2	1	-0.0001056284	-0.006306177	-0.0003942915	0.693	1
Hot_control.1_Acclimation vs Hot_control.6_Post-FMT2	1	0.0017235602	0.051851181	0.0032302306	0.746	1
Control.5_Post-FMT1 vs Treatment.5_Post-FMT1	1	0.1195408549	4.847647043	0.2442429086	0.063	1
Control.5_Post-FMT1 vs Hot_control.5_Post-FMT1	1	0.0525878365	1.773089316	0.0997625840	0.206	1
Control.5_Post-FMT1 vs Treatment.6_Post-FMT2	1	0.0265995825	1.175418056	0.0684360667	0.335	1
Control.5_Post-FMT1 vs Control.6_Post-FMT2	1	0.0145818992	0.699759916	0.0419023938	0.420	1
Control.5_Post-FMT1 vs Hot_control.6_Post-FMT2	1	-0.0080695208	-0.216173226	-0.0136958691	0.916	1
Treatment.5_Post-FMT1 vs Hot_control.5_Post-FMT1	1	0.0129803540	0.443076619	0.0286909552	0.460	1
Treatment.5_Post-FMT1 vs Treatment.6_Post-FMT2	1	0.0267162134	1.225605811	0.0755352882	0.311	1
Treatment.5_Post-FMT1 vs Control.6_Post-FMT2	1	0.0384388433	1.932815819	0.1141461550	0.220	1
Treatment.5_Post-FMT1 vs Hot_control.6_Post-FMT2	1	0.0553988290	1.478193905	0.0897060633	0.258	1
Hot_control.5_Post-FMT1 vs Treatment.6_Post-FMT2	1	-0.0040061386	-0.148504693	-0.0093684974	0.726	1
Hot_control.5_Post-FMT1 vs Control.6_Post-FMT2	1	0.0024023972	0.095389804	0.0059265296	0.612	1
Hot_control.5_Post-FMT1 vs Hot_control.6_Post-FMT2	1	-0.0004960759	-0.011903277	-0.0007445087	0.855	1
Treatment.6_Post-FMT2 vs Control.6_Post-FMT2	1	-0.0080428882	-0.442986255	-0.0284750185	0.859	1
Treatment.6_Post-FMT2 vs Hot_control.6_Post-FMT2	1	-0.0011796256	-0.034047378	-0.0021324990	0.897	1
Control.6_Post-FMT2 vs Hot_control.6_Post-FMT2	1	0.0036300838	0.110487573	0.0068581148	0.707	1

beta_richness_nmds_post6 <- richness_post6 %>%
  metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
  scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(subset_meta_post6, by = c("sample" = "Tube_code"))

beta_neutral_nmds_post6 <- neutral_post6 %>%
  metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
  scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(subset_meta_post6, by = c("sample" = "Tube_code"))

beta_phylo_nmds_post6 <- phylo_post6 %>%
  metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
  scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(subset_meta_post6, by = join_by(sample == Tube_code))

beta_func_nmds_post6 <- func_post6 %>%
  metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
  scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(subset_meta_post6, by = join_by(sample == Tube_code))

6.3.5 5. Are there differences between the control and the treatment group?

6.3.5.1 After 1 week –> Post-FMT1

6.3.5.1.1 Number of samples used

[1] 26

6.3.5.1.2 Richness

richness_post1 <- as.matrix(beta_q0n$S)
richness_post1 <- as.dist(richness_post1[rownames(richness_post1) %in% samples_to_keep_post1,
               colnames(richness_post1) %in% samples_to_keep_post1])
betadisper(richness_post1, subset_meta_post1$type) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     2 0.017675 0.0088373 2.3825    999  0.105
Residuals 23 0.085312 0.0037092                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
              Control Hot_control Treatment
Control                 0.0040000     0.668
Hot_control 0.0068795                 0.213
Treatment   0.6248469   0.2084296          

adonis2(richness_post1 ~ type,
        data = subset_meta_post1 %>% arrange(match(Tube_code,labels(richness_post1))),
        permutations = 999,
        strata = subset_meta_post1 %>% arrange(match(Tube_code,labels(richness_post1))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	2	1.195567	0.1448246	1.947534	1
Residual	23	7.059710	0.8551754	NA	NA
Total	25	8.255277	1.0000000	NA	NA

#Arrange of metadata dataframe
subset_meta_post1_arrange <- column_to_rownames(subset_meta_post1, "Tube_code")
subset_meta_post1_arrange<-subset_meta_post1_arrange[labels(richness_post1),]

pairwise <- pairwise.adonis(richness_post1, subset_meta_post1_arrange$type, perm=999)
pairwise%>%
  tt()

pairs	Df	SumsOfSqs	F.Model	R2	p.value	p.adjusted	sig
Control vs Treatment	1	0.5615418	1.729004	0.1033537	0.020	0.060
Control vs Hot_control	1	0.8438429	2.793772	0.1486541	0.001	0.003	*
Treatment vs Hot_control	1	0.3734921	1.268929	0.0779971	0.129	0.387

6.3.5.1.3 Neutral

neutral_post1 <- as.matrix(beta_q1n$S)
neutral_post1 <- as.dist(neutral_post1[rownames(neutral_post1) %in% samples_to_keep_post1,
               colnames(neutral_post1) %in% samples_to_keep_post1])
betadisper(neutral_post1, subset_meta_post1$type) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     2 0.011001 0.0055005 0.6303    999  0.572
Residuals 23 0.200714 0.0087267                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
            Control Hot_control Treatment
Control                 0.22500     0.955
Hot_control 0.21166                 0.446
Treatment   0.95468     0.43604          

adonis2(neutral_post1 ~ type,
        data = subset_meta_post1 %>% arrange(match(Tube_code,labels(neutral_post1))),
        permutations = 999,
        strata = subset_meta_post1 %>% arrange(match(Tube_code,labels(neutral_post1))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	2	1.395968	0.1900228	2.697931	1
Residual	23	5.950350	0.8099772	NA	NA
Total	25	7.346318	1.0000000	NA	NA

pairwise <- pairwise.adonis(neutral_post1, subset_meta_post1_arrange$type, perm=999)
pairwise%>%
  tt()

pairs	Df	SumsOfSqs	F.Model	R2	p.value	p.adjusted	sig
Control vs Treatment	1	0.6051778	2.250849	0.13047758	0.013	0.039	.
Control vs Hot_control	1	1.0528902	4.143637	0.20570451	0.002	0.006	*
Treatment vs Hot_control	1	0.4150076	1.637268	0.09840968	0.050	0.150

6.3.5.1.4 Phylogenetic

phylo_post1 <- as.matrix(beta_q1p$S)
phylo_post1 <- as.dist(phylo_post1[rownames(phylo_post1) %in% samples_to_keep_post1,
               colnames(phylo_post1) %in% samples_to_keep_post1])
betadisper(phylo_post1, subset_meta_post1$type) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     2 0.00440 0.0021994 0.1369    999  0.908
Residuals 23 0.36941 0.0160614                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
            Control Hot_control Treatment
Control                 0.93200     0.698
Hot_control 0.91505                 0.776
Treatment   0.63312     0.73046          

adonis2(phylo_post1 ~ type,
        data = subset_meta_post1 %>% arrange(match(Tube_code,labels(phylo_post1))),
        permutations = 999,
        strata = subset_meta_post1 %>% arrange(match(Tube_code,labels(phylo_post1))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	2	0.07451036	0.07059466	0.8735033	1
Residual	23	0.98095695	0.92940534	NA	NA
Total	25	1.05546731	1.00000000	NA	NA

pairwise <- pairwise.adonis(phylo_post1, subset_meta_post1_arrange$type, perm=999)
pairwise%>%
  tt()

pairs	Df	SumsOfSqs	F.Model	R2	p.value	p.adjusted	sig
Control vs Treatment	1	0.01842535	0.4144162	0.02688498	0.791	1.000
Control vs Hot_control	1	0.05987967	1.7387847	0.09802164	0.122	0.366
Treatment vs Hot_control	1	0.03212966	0.6477782	0.04139746	0.712	1.000

6.3.5.1.5 Functional

func_post1 <- as.matrix(beta_q1f$S)
func_post1 <- as.dist(func_post1[rownames(func_post1) %in% samples_to_keep_post1,
               colnames(func_post1) %in% samples_to_keep_post1])
betadisper(func_post1, subset_meta_post1$type) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     2 0.00400 0.0019999 0.1431    999   0.88
Residuals 23 0.32135 0.0139717                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
            Control Hot_control Treatment
Control                 0.61100     0.766
Hot_control 0.60188                 0.845
Treatment   0.74597     0.84473          

adonis2(func_post1 ~ type,
        data = subset_meta_post1 %>% arrange(match(Tube_code,labels(func_post1))),
        permutations = 999,
        strata = subset_meta_post1 %>% arrange(match(Tube_code,labels(func_post1))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	2	0.1230554	0.1608583	2.204479	1
Residual	23	0.6419374	0.8391417	NA	NA
Total	25	0.7649929	1.0000000	NA	NA

pairwise <- pairwise.adonis(func_post1, subset_meta_post1_arrange$type, perm=999)
pairwise%>%
  tt()

pairs	Df	SumsOfSqs	F.Model	R2	p.value	p.adjusted	sig
Control vs Treatment	1	0.11954085	4.8476470	0.24424291	0.086	0.258
Control vs Hot_control	1	0.05258784	1.7730893	0.09976258	0.245	0.735
Treatment vs Hot_control	1	0.01298035	0.4430766	0.02869096	0.452	1.000

beta_richness_nmds_post1 <- richness_post1 %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(subset_meta_post1, by = join_by(sample == Tube_code))

beta_neutral_nmds_post1 <- neutral_post1 %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(subset_meta_post1, by = join_by(sample == Tube_code))

beta_phylogenetic_nmds_post1 <- phylo_post1 %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(subset_meta_post1, by = join_by(sample == Tube_code))

beta_functional_nmds_post1 <- func_post1 %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(subset_meta_post1, by = join_by(sample == Tube_code))

p0<-beta_richness_nmds_post1 %>%
            group_by(type) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type)) +
        scale_color_manual(name="Type",
                       breaks=c("Control", "Hot_control", "Treatment"),
                       labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
                       values=c("#4477AA","#d57d2c","#76b183")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y = "NMDS2", x="NMDS1 \n Richness beta diversity") +
                theme_classic() +
                theme(legend.position="none")

p1<-beta_neutral_nmds_post1 %>%
            group_by(type) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type)) +
        scale_color_manual(name="Type",
                       breaks=c("Control", "Hot_control", "Treatment"),
                       labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
                       values=c("#4477AA","#d57d2c","#76b183")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y = "NMDS2", x="NMDS1 \n Neutral beta diversity") +
                theme_classic() +
                theme(legend.position="none")
  
p2<-beta_phylogenetic_nmds_post1 %>%
            group_by(type) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type)) +
                scale_color_manual(name="Type",
                       breaks=c("Control", "Hot_control", "Treatment"),
                       labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
                       values=c("#4477AA","#d57d2c","#76b183")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y= element_blank (), x="NMDS1 \n Phylogenetic beta diversity") +
                theme_classic() +
                theme(legend.position="none")

p3<-beta_functional_nmds_post1 %>%
            group_by(type) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type)) +
                scale_color_manual(name="Type",
                       breaks=c("Control", "Hot_control", "Treatment"),
                       labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
                       values=c("#4477AA","#d57d2c","#76b183")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y= element_blank (), x="NMDS1 \n Functional beta diversity") +
                theme_classic()+
                theme(legend.position="none")

ggarrange(p0, p1, p2, p3, ncol=2, nrow=2, common.legend = TRUE, legend="right")

6.3.5.2 After 2 weeks –>Post-FMT2

post2 <- meta %>%
  filter(time_point == "6_Post-FMT2")

post2.counts <- temp_genome_counts[,which(colnames(temp_genome_counts) %in% rownames(post2))]
identical(sort(colnames(post2.counts)),sort(as.character(rownames(post2))))

post2_nmds <- sample_metadata %>%
  filter(time_point == "6_Post-FMT2")

6.3.5.2.1 Number of samples used

[1] 27

6.3.5.2.2 Richness

richness_post2 <- as.matrix(beta_q0n$S)
richness_post2 <- as.dist(richness_post2[rownames(richness_post2) %in% samples_to_keep_post2,
               colnames(richness_post2) %in% samples_to_keep_post2])
betadisper(richness_post2, subset_meta_post2$type) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     2 0.002011 0.0010056 0.1982    999  0.826
Residuals 24 0.121775 0.0050740                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
            Control Hot_control Treatment
Control                 0.70400     0.801
Hot_control 0.67789                 0.645
Treatment   0.79246     0.59820          

adonis2(richness_post2 ~ type,
        data = subset_meta_post2 %>% arrange(match(Tube_code,labels(richness_post2))),
        permutations = 999,
        strata = subset_meta_post2 %>% arrange(match(Tube_code,labels(richness_post2))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	2	1.504341	0.1967776	2.939822	1
Residual	24	6.140538	0.8032224	NA	NA
Total	26	7.644879	1.0000000	NA	NA

#Arrange of metadata dataframe
subset_meta_post2_arrange <- column_to_rownames(subset_meta_post2, "Tube_code")
subset_meta_post2_arrange<-subset_meta_post2_arrange[labels(richness_post2),]

pairwise <- pairwise.adonis(richness_post2, subset_meta_post2_arrange$type, perm=999)
pairwise%>%
  tt()

pairs	Df	SumsOfSqs	F.Model	R2	p.value	p.adjusted	sig
Treatment vs Control	1	0.6463814	2.560441	0.1379515	0.002	0.006	*
Treatment vs Hot_control	1	0.4796256	1.916520	0.1069694	0.001	0.003	*
Control vs Hot_control	1	1.1305044	4.268317	0.2105906	0.001	0.003	*

6.3.5.2.3 Neutral

neutral_post2 <- as.matrix(beta_q1n$S)
neutral_post2 <- as.dist(neutral_post2[rownames(neutral_post2) %in% samples_to_keep_post2,
               colnames(neutral_post2) %in% samples_to_keep_post2])
betadisper(neutral_post2, subset_meta_post2$type) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     2 0.008262 0.0041311 0.8024    999  0.443
Residuals 24 0.123559 0.0051483                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
            Control Hot_control Treatment
Control                 0.44400     0.676
Hot_control 0.44675                 0.238
Treatment   0.65989     0.25095          

adonis2(neutral_post2 ~ type,
        data = subset_meta_post2 %>% arrange(match(Tube_code,labels(neutral_post2))),
        permutations = 999,
        strata = subset_meta_post2 %>% arrange(match(Tube_code,labels(neutral_post2))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	2	1.923807	0.2603795	4.224537	1
Residual	24	5.464666	0.7396205	NA	NA
Total	26	7.388473	1.0000000	NA	NA

pairwise <- pairwise.adonis(neutral_post2, subset_meta_post2_arrange$type, perm=999)
pairwise%>%
  tt()

pairs	Df	SumsOfSqs	F.Model	R2	p.value	p.adjusted	sig
Treatment vs Control	1	1.0227481	4.648335	0.2251191	0.001	0.003	*
Treatment vs Hot_control	1	0.5010202	2.206532	0.1211945	0.002	0.006	*
Control vs Hot_control	1	1.3619424	5.771031	0.2650785	0.001	0.003	*

6.3.5.2.4 Phylogenetic

phylo_post2 <- as.matrix(beta_q1p$S)
phylo_post2 <- as.dist(phylo_post2[rownames(phylo_post2) %in% samples_to_keep_post2,
               colnames(phylo_post2) %in% samples_to_keep_post2])
betadisper(phylo_post2, subset_meta_post2$type) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     2 0.000407 0.0002034 0.0487    999  0.961
Residuals 24 0.100305 0.0041794                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
            Control Hot_control Treatment
Control                 0.93700     0.852
Hot_control 0.93765                 0.757
Treatment   0.83933     0.76015          

adonis2(phylo_post2 ~ type,
        data = subset_meta_post2 %>% arrange(match(Tube_code,labels(phylo_post2))),
        permutations = 999,
        strata = subset_meta_post2 %>% arrange(match(Tube_code,labels(phylo_post2))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	2	0.1594363	0.2042241	3.079623	1
Residual	24	0.6212564	0.7957759	NA	NA
Total	26	0.7806927	1.0000000	NA	NA

pairwise <- pairwise.adonis(phylo_post2, subset_meta_post2_arrange$type, perm=999)
pairwise%>%
  tt()

pairs	Df	SumsOfSqs	F.Model	R2	p.value	p.adjusted	sig
Treatment vs Control	1	0.05927454	2.382025	0.1295845	0.032	0.096
Treatment vs Hot_control	1	0.06906280	2.722460	0.1454115	0.005	0.015	.
Control vs Hot_control	1	0.11081709	4.043656	0.2017424	0.001	0.003	*

6.3.5.2.5 Functional

func_post2 <- as.matrix(beta_q1f$S)
func_post2 <- as.dist(func_post2[rownames(func_post2) %in% samples_to_keep_post2,
               colnames(func_post2) %in% samples_to_keep_post2])
betadisper(func_post2, subset_meta_post2$type) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     2 0.01259 0.0062962 0.3249    999   0.79
Residuals 24 0.46507 0.0193778                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
            Control Hot_control Treatment
Control                 0.53000     0.638
Hot_control 0.45381                 0.774
Treatment   0.57452     0.74365          

adonis2(func_post2 ~ type,
        data = subset_meta_post2 %>% arrange(match(Tube_code,labels(func_post2))),
        permutations = 999,
        strata = subset_meta_post2 %>% arrange(match(Tube_code,labels(func_post2))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	2	-0.003728287	-0.005470434	-0.06528805	1
Residual	24	0.685262325	1.005470434	NA	NA
Total	26	0.681534039	1.000000000	NA	NA

pairwise <- pairwise.adonis(func_post2, subset_meta_post2_arrange$type, perm=999)
pairwise%>%
  tt()

pairs	Df	SumsOfSqs	F.Model	R2	p.value	p.adjusted	sig
Treatment vs Control	1	-0.008042888	-0.44298625	-0.028475019	0.829	1
Treatment vs Hot_control	1	-0.001179626	-0.03404738	-0.002132499	0.890	1
Control vs Hot_control	1	0.003630084	0.11048757	0.006858115	0.693	1

beta_richness_nmds_post2 <- richness_post2 %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(subset_meta_post2, by = join_by(sample == Tube_code))

beta_neutral_nmds_post2 <- neutral_post2 %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(subset_meta_post2, by = join_by(sample == Tube_code))

beta_phylogenetic_nmds_post2 <- phylo_post2 %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(subset_meta_post2, by = join_by(sample == Tube_code))

beta_functional_nmds_post2 <- func_post2 %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(subset_meta_post2, by = join_by(sample == Tube_code))

p0<-beta_richness_nmds_post2 %>%
            group_by(type) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type)) +
        scale_color_manual(name="Type",
                       breaks=c("Control", "Hot_control", "Treatment"),
                       labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
                       values=c("#4477AA","#d57d2c","#76b183")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y = "NMDS2", x="NMDS1 \n Richness beta diversity") +
                theme_classic() +
                theme(legend.position="none")

p1<-beta_neutral_nmds_post2 %>%
            group_by(type) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type)) +
        scale_color_manual(name="Type",
                       breaks=c("Control", "Hot_control", "Treatment"),
                       labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
                       values=c("#4477AA","#d57d2c","#76b183")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y = "NMDS2", x="NMDS1 \n Neutral beta diversity") +
                theme_classic() +
                theme(legend.position="none")
  
p2<-beta_phylogenetic_nmds_post2 %>%
            group_by(type) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type)) +
                scale_color_manual(name="Type",
                       breaks=c("Control", "Hot_control", "Treatment"),
                       labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
                       values=c("#4477AA","#d57d2c","#76b183")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y= element_blank (), x="NMDS1 \n Phylogenetic beta diversity") +
                theme_classic() +
                theme(legend.position="none")

p3<-beta_functional_nmds_post2 %>%
            group_by(type) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type)) +
                scale_color_manual(name="Type",
                       breaks=c("Control", "Hot_control", "Treatment"),
                       labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
                       values=c("#4477AA","#d57d2c","#76b183")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y= element_blank (), x="NMDS1 \n Functional beta diversity") +
                theme_classic()+
                theme(legend.position="none")

ggarrange(p0, p1, p2, p3, ncol=2, nrow=2, common.legend = TRUE, legend="right")

6.3.5.3 Post1 vs Post2

post5 <- meta %>%
  filter(time_point == "6_Post-FMT2" | time_point == "5_Post-FMT1")

post5.counts <- temp_genome_counts[,which(colnames(temp_genome_counts) %in% rownames(post5))]
identical(sort(colnames(post5.counts)),sort(as.character(rownames(post5))))

post5_nmds <- sample_metadata %>%
  filter(time_point == "6_Post-FMT2"| time_point == "5_Post-FMT1")

6.3.5.3.1 Number of samples used

[1] 53

6.3.5.3.2 Richness

richness_post5 <- as.matrix(beta_q0n$S)
richness_post5 <- as.dist(richness_post5[rownames(richness_post5) %in% samples_to_keep_post5,
               colnames(richness_post5) %in% samples_to_keep_post5])
betadisper(richness_post5, subset_meta_post5$type) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     2 0.01841 0.0092048 1.7364    999   0.19
Residuals 50 0.26505 0.0053010                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
             Control Hot_control Treatment
Control                 0.037000     0.712
Hot_control 0.039117                 0.227
Treatment   0.716358    0.218648          

adonis2(richness_post5 ~ type*time_point,
        data = subset_meta_post5 %>% arrange(match(Tube_code,labels(richness_post5))),
        permutations = 999,
        strata = subset_meta_post5 %>% arrange(match(Tube_code,labels(richness_post5))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	5	3.261916	0.1981463	2.322836	0.001
Residual	47	13.200248	0.8018537	NA	NA
Total	52	16.462164	1.0000000	NA	NA

#Arrange of metadata dataframe
subset_meta_post5_arrange <- column_to_rownames(subset_meta_post5, "Tube_code")
subset_meta_post5_arrange<-subset_meta_post5_arrange[labels(richness_post5),]

pairwise <- pairwise.adonis(richness_post5, subset_meta_post5_arrange$type_time, perm=999)
pairwise%>%
  tt()

pairs	Df	SumsOfSqs	F.Model	R2	p.value	p.adjusted	sig
Control.5_Post-FMT1 vs Treatment.5_Post-FMT1	1	0.5615418	1.729004	0.10335366	0.015	0.225
Control.5_Post-FMT1 vs Hot_control.5_Post-FMT1	1	0.8438429	2.793772	0.14865413	0.001	0.015	.
Control.5_Post-FMT1 vs Treatment.6_Post-FMT2	1	0.7628135	2.683925	0.14364890	0.001	0.015	.
Control.5_Post-FMT1 vs Control.6_Post-FMT2	1	0.3432605	1.148733	0.06698647	0.264	1.000
Control.5_Post-FMT1 vs Hot_control.6_Post-FMT2	1	1.1269580	3.799256	0.19188884	0.001	0.015	.
Treatment.5_Post-FMT1 vs Hot_control.5_Post-FMT1	1	0.3734921	1.268929	0.07799710	0.097	1.000
Treatment.5_Post-FMT1 vs Treatment.6_Post-FMT2	1	0.3571397	1.297184	0.07959561	0.137	1.000
Treatment.5_Post-FMT1 vs Control.6_Post-FMT2	1	0.7769467	2.670898	0.15114670	0.001	0.015	.
Treatment.5_Post-FMT1 vs Hot_control.6_Post-FMT2	1	0.6502360	2.253407	0.13060650	0.002	0.030	.
Hot_control.5_Post-FMT1 vs Treatment.6_Post-FMT2	1	0.4132091	1.616138	0.09174188	0.012	0.180
Hot_control.5_Post-FMT1 vs Control.6_Post-FMT2	1	1.0163992	3.760571	0.19030682	0.001	0.015	.
Hot_control.5_Post-FMT1 vs Hot_control.6_Post-FMT2	1	0.2732563	1.019281	0.05988979	0.433	1.000
Treatment.6_Post-FMT2 vs Control.6_Post-FMT2	1	0.6463814	2.560441	0.13795154	0.001	0.015	.
Treatment.6_Post-FMT2 vs Hot_control.6_Post-FMT2	1	0.4796256	1.916520	0.10696943	0.001	0.015	.
Control.6_Post-FMT2 vs Hot_control.6_Post-FMT2	1	1.1305044	4.268317	0.21059061	0.001	0.015	.

6.3.5.3.3 Neutral

neutral_post5 <- as.matrix(beta_q1n$S)
neutral_post5 <- as.dist(neutral_post5[rownames(neutral_post5) %in% samples_to_keep_post5,
               colnames(neutral_post5) %in% samples_to_keep_post5])
betadisper(neutral_post5, subset_meta_post5$type) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq   Mean Sq     F N.Perm Pr(>F)
Groups     2 0.01992 0.0099587 1.565    999  0.212
Residuals 50 0.31818 0.0063636                    

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
            Control Hot_control Treatment
Control                 0.11400     0.876
Hot_control 0.10701                 0.165
Treatment   0.87156     0.17449          

adonis2(neutral_post5 ~ type*time_point,
        data = subset_meta_post5 %>% arrange(match(Tube_code,labels(neutral_post5))),
        permutations = 999,
        strata = subset_meta_post5 %>% arrange(match(Tube_code,labels(neutral_post5))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	5	3.947979	0.2569798	3.251069	0.001
Residual	47	11.415016	0.7430202	NA	NA
Total	52	15.362995	1.0000000	NA	NA

pairwise <- pairwise.adonis(neutral_post5, subset_meta_post5_arrange$type_time, perm=999)
pairwise%>%
  tt()

pairs	Df	SumsOfSqs	F.Model	R2	p.value	p.adjusted	sig
Control.5_Post-FMT1 vs Treatment.5_Post-FMT1	1	0.6051778	2.250849	0.13047758	0.019	0.285
Control.5_Post-FMT1 vs Hot_control.5_Post-FMT1	1	1.0528902	4.143637	0.20570451	0.001	0.015	.
Control.5_Post-FMT1 vs Treatment.6_Post-FMT2	1	0.8908158	3.714692	0.18842252	0.002	0.030	.
Control.5_Post-FMT1 vs Control.6_Post-FMT2	1	0.3860927	1.552176	0.08843210	0.094	1.000
Control.5_Post-FMT1 vs Hot_control.6_Post-FMT2	1	1.3122237	5.130273	0.24279254	0.002	0.030	.
Treatment.5_Post-FMT1 vs Hot_control.5_Post-FMT1	1	0.4150076	1.637268	0.09840968	0.056	0.840
Treatment.5_Post-FMT1 vs Treatment.6_Post-FMT2	1	0.3157079	1.325203	0.08117526	0.152	1.000
Treatment.5_Post-FMT1 vs Control.6_Post-FMT2	1	1.0579520	4.270010	0.22158835	0.002	0.030	.
Treatment.5_Post-FMT1 vs Hot_control.6_Post-FMT2	1	0.7454015	2.920049	0.16294873	0.001	0.015	.
Hot_control.5_Post-FMT1 vs Treatment.6_Post-FMT2	1	0.4377161	1.942126	0.10824392	0.006	0.090
Hot_control.5_Post-FMT1 vs Control.6_Post-FMT2	1	1.3766597	5.875279	0.26858075	0.001	0.015	.
Hot_control.5_Post-FMT1 vs Hot_control.6_Post-FMT2	1	0.3176516	1.316137	0.07600637	0.199	1.000
Treatment.6_Post-FMT2 vs Control.6_Post-FMT2	1	1.0227481	4.648335	0.22511910	0.002	0.030	.
Treatment.6_Post-FMT2 vs Hot_control.6_Post-FMT2	1	0.5010202	2.206532	0.12119453	0.001	0.015	.
Control.6_Post-FMT2 vs Hot_control.6_Post-FMT2	1	1.3619424	5.771031	0.26507845	0.001	0.015	.

6.3.5.3.4 Phylogenetic

phylo_post5 <- as.matrix(beta_q1p$S)
phylo_post5 <- as.dist(phylo_post5[rownames(phylo_post5) %in% samples_to_keep_post5,
               colnames(phylo_post5) %in% samples_to_keep_post5])
betadisper(phylo_post5, subset_meta_post5$type) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     2 0.00051 0.0002543 0.0265    999  0.973
Residuals 50 0.47996 0.0095993                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
            Control Hot_control Treatment
Control                 0.88500     0.837
Hot_control 0.88926                 0.926
Treatment   0.82391     0.91902          

adonis2(phylo_post5 ~ type*time_point,
        data = subset_meta_post5 %>% arrange(match(Tube_code,labels(phylo_post5))),
        permutations = 999,
        strata = subset_meta_post5 %>% arrange(match(Tube_code,labels(phylo_post5))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	5	0.3518222	0.180049	2.0641	0.001
Residual	47	1.6022134	0.819951	NA	NA
Total	52	1.9540356	1.000000	NA	NA

pairwise <- pairwise.adonis(phylo_post5, subset_meta_post5_arrange$type_time, perm=999)
pairwise%>%
  tt()

pairs	Df	SumsOfSqs	F.Model	R2	p.value	p.adjusted	sig
Control.5_Post-FMT1 vs Treatment.5_Post-FMT1	1	0.01842535	0.4144162	0.02688498	0.768	1.000
Control.5_Post-FMT1 vs Hot_control.5_Post-FMT1	1	0.05987967	1.7387847	0.09802164	0.134	1.000
Control.5_Post-FMT1 vs Treatment.6_Post-FMT2	1	0.07917244	3.0180046	0.15869197	0.010	0.150
Control.5_Post-FMT1 vs Control.6_Post-FMT2	1	0.04335491	1.5335604	0.08746429	0.196	1.000
Control.5_Post-FMT1 vs Hot_control.6_Post-FMT2	1	0.10783045	3.7500438	0.18987521	0.002	0.030	.
Treatment.5_Post-FMT1 vs Hot_control.5_Post-FMT1	1	0.03212966	0.6477782	0.04139746	0.689	1.000
Treatment.5_Post-FMT1 vs Treatment.6_Post-FMT2	1	0.06393539	1.5651817	0.09448624	0.131	1.000
Treatment.5_Post-FMT1 vs Control.6_Post-FMT2	1	0.05265949	1.2240203	0.07544494	0.295	1.000
Treatment.5_Post-FMT1 vs Hot_control.6_Post-FMT2	1	0.09753501	2.2402429	0.12994265	0.016	0.240
Hot_control.5_Post-FMT1 vs Treatment.6_Post-FMT2	1	0.07228545	2.3279593	0.12701683	0.031	0.465
Hot_control.5_Post-FMT1 vs Control.6_Post-FMT2	1	0.11759094	3.5538444	0.18174658	0.001	0.015	.
Hot_control.5_Post-FMT1 vs Hot_control.6_Post-FMT2	1	0.06667255	1.9859527	0.11041687	0.097	1.000
Treatment.6_Post-FMT2 vs Control.6_Post-FMT2	1	0.05927454	2.3820253	0.12958449	0.025	0.375
Treatment.6_Post-FMT2 vs Hot_control.6_Post-FMT2	1	0.06906280	2.7224602	0.14541146	0.003	0.045	.
Control.6_Post-FMT2 vs Hot_control.6_Post-FMT2	1	0.11081709	4.0436561	0.20174244	0.002	0.030	.

6.3.5.3.5 Functional

func_post5 <- as.matrix(beta_q1f$S)
func_post5 <- as.dist(func_post5[rownames(func_post5) %in% samples_to_keep_post5,
               colnames(func_post5) %in% samples_to_keep_post5])
betadisper(func_post5, subset_meta_post5$type) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     2 0.00785 0.0039232 0.2322    999  0.818
Residuals 50 0.84483 0.0168966                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
            Control Hot_control Treatment
Control                 0.57000     0.599
Hot_control 0.52384                 0.860
Treatment   0.58787     0.85068          

adonis2(func_post5 ~ type*time_point,
        data = subset_meta_post5 %>% arrange(match(Tube_code,labels(func_post5))),
        permutations = 999,
        strata = subset_meta_post5 %>% arrange(match(Tube_code,labels(func_post5))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
Model	5	0.1076682	0.07503698	0.7625685	0.499
Residual	47	1.3271997	0.92496302	NA	NA
Total	52	1.4348679	1.00000000	NA	NA

pairwise <- pairwise.adonis(func_post5, subset_meta_post5_arrange$type_time, perm=999)
pairwise%>%
  tt()

pairs	Df	SumsOfSqs	F.Model	R2	p.value	p.adjusted	sig
Control.5_Post-FMT1 vs Treatment.5_Post-FMT1	1	0.1195408549	4.84764704	0.2442429086	0.072	1
Control.5_Post-FMT1 vs Hot_control.5_Post-FMT1	1	0.0525878365	1.77308932	0.0997625840	0.218	1
Control.5_Post-FMT1 vs Treatment.6_Post-FMT2	1	0.0265995825	1.17541806	0.0684360667	0.303	1
Control.5_Post-FMT1 vs Control.6_Post-FMT2	1	0.0145818992	0.69975992	0.0419023938	0.437	1
Control.5_Post-FMT1 vs Hot_control.6_Post-FMT2	1	-0.0080695208	-0.21617323	-0.0136958691	0.921	1
Treatment.5_Post-FMT1 vs Hot_control.5_Post-FMT1	1	0.0129803540	0.44307662	0.0286909552	0.480	1
Treatment.5_Post-FMT1 vs Treatment.6_Post-FMT2	1	0.0267162134	1.22560581	0.0755352882	0.338	1
Treatment.5_Post-FMT1 vs Control.6_Post-FMT2	1	0.0384388433	1.93281582	0.1141461550	0.239	1
Treatment.5_Post-FMT1 vs Hot_control.6_Post-FMT2	1	0.0553988290	1.47819391	0.0897060633	0.260	1
Hot_control.5_Post-FMT1 vs Treatment.6_Post-FMT2	1	-0.0040061386	-0.14850469	-0.0093684974	0.735	1
Hot_control.5_Post-FMT1 vs Control.6_Post-FMT2	1	0.0024023972	0.09538980	0.0059265296	0.604	1
Hot_control.5_Post-FMT1 vs Hot_control.6_Post-FMT2	1	-0.0004960759	-0.01190328	-0.0007445087	0.849	1
Treatment.6_Post-FMT2 vs Control.6_Post-FMT2	1	-0.0080428882	-0.44298625	-0.0284750185	0.832	1
Treatment.6_Post-FMT2 vs Hot_control.6_Post-FMT2	1	-0.0011796256	-0.03404738	-0.0021324990	0.913	1
Control.6_Post-FMT2 vs Hot_control.6_Post-FMT2	1	0.0036300838	0.11048757	0.0068581148	0.703	1

beta_richness_nmds_post5 <- richness_post5 %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(subset_meta_post5, by = join_by(sample == Tube_code))

beta_neutral_nmds_post5 <- neutral_post5 %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(subset_meta_post5, by = join_by(sample == Tube_code))

beta_phylogenetic_nmds_post5 <- phylo_post5 %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(subset_meta_post5, by = join_by(sample == Tube_code))

beta_functional_nmds_post5 <- func_post5 %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(subset_meta_post5, by = join_by(sample == Tube_code))

p0<-beta_richness_nmds_post5 %>%
            group_by(type, time_point) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
        scale_color_manual(name="Type",
                       breaks=c("Control", "Hot_control", "Treatment"),
                       labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
                       values=c("#4477AA","#d57d2c","#76b183")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y = "NMDS2", x="NMDS1 \n Richness beta diversity") +
                theme_classic() +
                theme(legend.position="none")

p1<-beta_neutral_nmds_post5 %>%
            group_by(type, time_point) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
        scale_color_manual(name="Type",
                       breaks=c("Control", "Hot_control", "Treatment"),
                       labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
                       values=c("#4477AA","#d57d2c","#76b183")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y = "NMDS2", x="NMDS1 \n Neutral beta diversity") +
                theme_classic() +
                theme(legend.position="none")
  
p2<-beta_phylogenetic_nmds_post5 %>%
            group_by(type, time_point) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
                scale_color_manual(name="Type",
                       breaks=c("Control", "Hot_control", "Treatment"),
                       labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
                       values=c("#4477AA","#d57d2c","#76b183")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y= element_blank (), x="NMDS1 \n Phylogenetic beta diversity") +
                theme_classic() +
                theme(legend.position="none")

p3<-beta_functional_nmds_post5 %>%
            group_by(type, time_point) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
                scale_color_manual(name="Type",
                       breaks=c("Control", "Hot_control", "Treatment"),
                       labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
                       values=c("#4477AA","#d57d2c","#76b183")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y= element_blank (), x="NMDS1 \n Functional beta diversity") +
                theme_classic()+
                theme(legend.position="none")

ggarrange(p0, p1, p2, p3, ncol=2, nrow=2, common.legend = TRUE, legend="right")